Libraries

Include

using Revise

Algorithms

using DifferentialEquations
using Optim
using Statistics
using Distances

Data

using CSV
using DataFrames
using DataFramesMeta
using Chain: @chain
using Pipe: @pipe
using DataStructures

Plotting

using Plots
using StatsPlots
pyplot()

Plots.PyPlotBackend()

Modelling framework

includet("models/base.jl")

Data

We don't need much data, but having something to put in our plots would be handy.

Experimental data

data = map(["wt", "ifnko"]) do geno
    active = @linq CSV.read("data_old/data_$(geno)_active.csv", DataFrame) |>
        transform(name="active") |>
        select(genotype=:genotype, age=:age, name=:name, value=:active)
    counts = @linq CSV.read("data_old/data_$(geno)_counts.csv", DataFrame) |>
        transform(name="total") |>
        select(genotype=:genotype, age=:age, name=:name, value=:count)
    vcat([active, counts]...)
end
data = vcat(data...)

┌ Warning: Using an un-quoted Symbol on the LHS is deprecated. Write :name = ... instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/parsing.jl:298
┌ Warning: Using an un-quoted Symbol on the LHS is deprecated. Write :genotype = ... instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/parsing.jl:298
┌ Warning: Using an un-quoted Symbol on the LHS is deprecated. Write :age = ... instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/parsing.jl:298
┌ Warning: Using an un-quoted Symbol on the LHS is deprecated. Write :name = ... instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/parsing.jl:298
┌ Warning: Using an un-quoted Symbol on the LHS is deprecated. Write :value = ... instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/parsing.jl:298
┌ Warning: Using an un-quoted Symbol on the LHS is deprecated. Write :name = ... instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/parsing.jl:298
┌ Warning: Using an un-quoted Symbol on the LHS is deprecated. Write :genotype = ... instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/parsing.jl:298
┌ Warning: Using an un-quoted Symbol on the LHS is deprecated. Write :age = ... instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/parsing.jl:298
┌ Warning: Using an un-quoted Symbol on the LHS is deprecated. Write :name = ... instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/parsing.jl:298
┌ Warning: Using an un-quoted Symbol on the LHS is deprecated. Write :value = ... instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/parsing.jl:298

41 rows × 4 columns

	genotype	age	name	value
	String7	Int64	String	Float64
1	wt	660	active	0.234043
2	wt	660	active	0.15
3	wt	660	active	0.2
4	wt	660	active	0.228571
5	wt	60	active	0.586207
6	wt	60	active	0.467742
7	wt	60	active	0.280702
8	wt	60	total	764.0
9	wt	60	total	755.5
10	wt	60	total	1254.5
11	wt	60	total	845.5
12	wt	60	total	708.0
13	wt	60	total	1221.0
14	wt	60	total	804.0
15	wt	60	total	699.0
16	wt	210	total	232.0
17	wt	210	total	396.0
18	wt	210	total	238.0
19	wt	660	total	116.333
20	wt	660	total	88.3333
21	wt	660	total	226.0
22	ifnko	660	active	0.392857
23	ifnko	660	active	0.428571
24	ifnko	660	active	0.35
25	ifnko	60	active	0.333333
26	ifnko	60	active	0.352941
27	ifnko	60	active	0.340909
28	ifnko	60	total	584.0
29	ifnko	60	total	691.0
30	ifnko	60	total	1022.5
⋮	⋮	⋮	⋮	⋮

Simulation data

sim_ifnko_both = CSV.read("data_old/solution_ifnko_self_renewal_and_quiescence.csv", DataFrame)
sim_wt_both = CSV.read("data_old/solution_wt_self_renewal_and_quiescence.csv", DataFrame)
sim = vcat(@transform(sim_wt_both, genotype="wt"), @transform(sim_ifnko_both, genotype="ko"))

┌ Warning: Using an un-quoted Symbol on the LHS is deprecated. Write :genotype = ... instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/parsing.jl:298
┌ Warning: Using an un-quoted Symbol on the LHS is deprecated. Write :genotype = ... instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/parsing.jl:298

1,087 rows × 6 columns

	t	active	counts	b	r	genotype
	Float64	Float64	Float64	Float64	Float64	String
1	0.0	0.474127	1986.31	0.482161	0.842565	wt
2	0.157915	0.474119	1981.27	0.482174	0.842313	wt
3	0.271147	0.474101	1977.67	0.482184	0.842133	wt
4	0.423967	0.474068	1972.82	0.482197	0.84189	wt
5	0.599921	0.474017	1967.26	0.482212	0.84161	wt
6	0.804687	0.473946	1960.82	0.48223	0.841284	wt
7	1.0378	0.473854	1953.52	0.48225	0.840913	wt
8	1.29979	0.473741	1945.36	0.482273	0.840497	wt
9	1.59052	0.473609	1936.35	0.482298	0.840035	wt
10	1.90964	0.473459	1926.54	0.482325	0.839528	wt
11	2.25654	0.473291	1915.94	0.482355	0.838978	wt
12	2.63039	0.473109	1904.62	0.482387	0.838385	wt
13	3.03023	0.472912	1892.6	0.482421	0.837752	wt
14	3.45497	0.472702	1879.95	0.482458	0.837079	wt
15	3.90343	0.47248	1866.72	0.482496	0.83637	wt
16	4.37437	0.472247	1852.97	0.482536	0.835626	wt
17	4.86656	0.472002	1838.74	0.482578	0.834848	wt
18	5.37872	0.471748	1824.09	0.482621	0.83404	wt
19	5.9096	0.471485	1809.08	0.482666	0.833204	wt
20	6.458	0.471213	1793.75	0.482712	0.83234	wt
21	7.02271	0.470933	1778.15	0.482759	0.831452	wt
22	7.6026	0.470646	1762.33	0.482808	0.830541	wt
23	8.19658	0.470351	1746.32	0.482857	0.829609	wt
24	8.80361	0.470051	1730.17	0.482908	0.828658	wt
25	9.42271	0.469744	1713.91	0.482959	0.827688	wt
26	10.053	0.469432	1697.58	0.483011	0.826703	wt
27	10.6935	0.469115	1681.19	0.483064	0.825702	wt
28	11.3435	0.468793	1664.8	0.483118	0.824688	wt
29	12.0023	0.468467	1648.4	0.483172	0.823662	wt
30	12.669	0.468138	1632.04	0.483226	0.822624	wt
⋮	⋮	⋮	⋮	⋮	⋮	⋮

Motivation plot

This plot visualises the problem. Due to the specific dynamics in the niche wt has more output from stem cells early in life. Past 400 days of age however, due to the lack of increased quiescence, knockout mice have much more active stem cells and thus higher output.

p = plot(yscale=:log10, ylab="rate of TAP creation by NSC division", xlab="Time (days)")
for group in groupby(sim, :genotype)
    group = @transform(group, tapo=(1 .- :b) .* :counts .* :active .* 2 * 0.95)
    @df group plot!(p, :t, :tapo, label=unique(group.genotype)[1], lw=2)
end
p

┌ Warning: Using an un-quoted Symbol on the LHS is deprecated. Write :tapo = ... instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/parsing.jl:298

@chain sim begin
    @transform tapo=(1 .- :b) .* :counts .* :active .* 2 * 0.95 
    groupby(:genotype)
    @combine mean=mean(:tapo)
end

┌ Warning: Using an un-quoted Symbol on the LHS is deprecated. Write :tapo = ... instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/parsing.jl:298
┌ Warning: Using an un-quoted Symbol on the LHS is deprecated. Write :mean = ... instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/parsing.jl:298

2 rows × 2 columns

	genotype	mean
	String	Float64
1	wt	141.412
2	ko	228.919

Regular Kalamakis model

Definition

includet("models/nsc_kalamakis.jl")

times = collect(1:0.1:700.0)
model = KalamakisFull((0.0, 700.0), times, zeros(Int64, length(times)), Dict(:pₛ => 0.95))

KalamakisFull((0.0, 700.0), [1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9  …  699.1, 699.2, 699.3, 699.4, 699.5, 699.6, 699.7, 699.8, 699.9, 700.0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0  …  0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [6 => 0.95])

parameter_names(model)

6-element Vector{Symbol}:
 :nsc₀
 :r₀
 :βr
 :b₀
 :βb
 :pₛ

Previously estimated parameters

param = Dict(
    :wt => Dict(
        :nsc₀ => 1986.0,
        :r₀ => 0.8426,
        :βr => 0.001891,
        :b₀ => 0.4822,
        :βb => 0.004862
    ),
    :ko => Dict(
        :nsc₀ => 2500.0,
        :r₀ => 0.4928,
        :βr => 1e-12,
        :b₀ => 0.4455,
        :βb => 0.01288
    )
)

Dict{Symbol, Dict{Symbol, Float64}} with 2 entries:
  :wt => Dict(:βb=>0.004862, :b₀=>0.4822, :βr=>0.001891, :r₀=>0.8426, :nsc₀=>19…
  :ko => Dict(:βb=>0.01288, :b₀=>0.4455, :βr=>1.0e-12, :r₀=>0.4928, :nsc₀=>2500…

Solutions

sim = simulate(model, param[:wt]) |> link(model)
ylabs = link_names(model)
p1 = plot(times, sim[:,1], lab="", lc=:red, lw=2, title="wt", yscale=:log10, ylab=ylabs[1])
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
p2 = plot(times, sim[:,2], lab="", lc=:red, lw=2, ylab=ylabs[2])
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
plot(p1, p2, xlab="Time (days)")

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

sim = simulate(model, param[:ko]) |> link(model)
ylabs = link_names(model)
p1 = plot(times, sim[:,1], lab="", lc=:red, lw=2, title="ko", yscale=:log10, ylab=ylabs[1])
@df @where(data, :genotype .== "ifnko", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
p2 = plot(times, sim[:,2], lab="", lc=:red, lw=2, ylab=ylabs[2])
@df @where(data, :genotype .== "ifnko", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
plot(p1, p2, xlab="Time (days)")

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

┌ Warning: No strict ticks found
└ @ PlotUtils /home/jooa/.julia/packages/PlotUtils/VgXdq/src/ticks.jl:294

Parameters

param_rt = Dict(
    geno => exponential_decay(params[:r₀], params[:βr])
    for (geno, params) in param
)
param_bt = Dict(
    geno => flattening_curve(params[:b₀], params[:βb])
    for (geno, params) in param
)

Dict{Symbol, var"#64#65"{Float64, Float64}} with 2 entries:
  :wt => #64
  :ko => #64

p1 = plot(times, param_rt[:wt], lab="wt", ylab="Activation rate")
plot!(p1, times, param_rt[:ko], lab="ko")
p2 = plot(times, param_bt[:wt], lab="wt", ylab="Self-renewal probability", legend=:bottomright)
plot!(p2, times, param_bt[:ko], lab="ko")
plot(p1, p2, xlab="Time (days)", lw=2)

Differences from other notebook.

We take the intervention models that we previously built. Further, we assume that $b$ is a function $\text{NSC}\rightarrow b$ with two parameters, one of which is dependent on genotype. More specifically, $$ b(\text{NSC}) = \frac{1}{1+\left(\frac{K_{ab}}{\text{NSC}}\right)^{n_b}} $$ Where $K_{ab}$ is the shared half-probability population around 150 stem cells and $n_b$ the genotype dependent shape parameter.

hill(ka, n, x) = 1/(1+(ka/x)^n)
hill(ka, n) = x -> hill(ka, n, x)

hill (generic function with 2 methods)

A place to store metrics

metrics_neurogen = DataFrame(itime = collect(0:1:660), type="neurogenesis");
metrics_loss = DataFrame(itime = collect(0:1:660), type="aged loss");

Beta/self model

Takes parameters, outputs simulations. Plots simulations with datapoints. Thus needs:

• Data
• Parameters

In this model the effect of interferon is a changed $r'$, as well as a changedl $b_n$. This is assuming that in wt: $$ r'_\text{wt}(t) \neq 0, $$ i.e. $ \beta_r \neq 0. $ and in ko $$ r'_\text{ko}(t) = 0, $$ i.e. $ \beta_r = 0. $

Since $r$ is exponential this means that this model should for some intervention timepoint $t_1$ have a $r$ of the form: $$ \frac{\text{d}}{\text{d}t}r(t) = \begin{cases} -\beta_r r(t)& \text{if }t<t_1 \\ 0 & \text{else} \end{cases} $$

Thus we have parameters: $$ \theta = ( \text{NSC}_0, r_0, \beta_r, K_{ab}, n_{b,wt}, n_{b,ko}, p_s, t_1 ) $$

struct ChangedBetaSelf <: ODEModel
    tspan::Tuple{Float64, Float64}
    times::Vector{Float64}
    values::Vector{Int64}
    fixed::Vector{Pair{Int64, Float64}}
    ChangedBetaSelf(tspan, times, values, fixed::Vector{Pair{I, N}}) where {I <: Integer, N <: Number} = new(tspan, times, values, fixed)
end
function ratefun(model::ChangedBetaSelf)
    function(du, u, p, t)
        _, _, βr, kab_wt, kab_ko, nb_wt, nb_ko, pₛ, t₁ = p
        Q = u[1]; A = u[2]; r = u[3]
        b = 1 - hill(t < t₁ ? kab_wt : kab_ko, 
                     t < t₁ ? nb_wt : nb_ko, 
                     Q+A)
        du[1] = dQ = -r * Q + 2* b * pₛ * A
        du[2] = dA = r * Q - pₛ * A
        if t < t₁
            du[3] = dr = -βr * r
        else
            du[3] = dr = 0
        end
    end
end
function initial(t::ChangedBetaSelf, x::AbstractVector)
    nsc₀, r₀, _, kab_wt, _, nb_wt, _, pₛ, _ = x
    b₀ = hill(kab_wt, nb_wt, 0)
    ratio = sqrt(((pₛ - r₀)/(2*r₀))^2 + (2*b₀*pₛ) / r₀)
    vcat(nsc₀ .* [1-1/(ratio+1), 1/(ratio+1)], [r₀])
end
link(t::ChangedBetaSelf, x::AbstractArray) = hcat(x[:,1] .+ x[:,2], x[:,2] ./ (x[:,1] .+ x[:,2]))
parameter_names(t::Type{ChangedBetaSelf}) = [:nsc₀, :r₀, :βr, :kab_wt, :kab_ko, :nb_wt, :nb_ko, :pₛ, :t₁]
bounds(t::ChangedBetaSelf) = [(100.0, 10000.0), (0.0, 1.0), (0.0, 0.1), (0.0, 1000.0), (0.0, 1000.0), (0.01, 0.1), (0.01, 0.1), (0.0, 1.0), (0.0, 700.0)]
output_names(t::ChangedBetaSelf) = ["qNSC", "aNSC", "r"]
link_names(t::ChangedBetaSelf) = ["Total NSC", "Fraction active NSC"]

link_names (generic function with 3 methods)

Self-renewal depends on genotype

param_beta = Dict( # from partial cofit
    #:nsc₀ => 1986.0,
    :nsc₀ => 3000.0,
    #:r₀ => 0.8426,
    :r₀ => 0.806991,
    #:βr => 0.001891,
    :βr => 0.00184341,
    :kab_wt => 155.445,
    :kab_ko => 148.364,
    :nb_wt => 0.0399102,
    :nb_ko => 0.0843197,
    :pₛ => 0.95,
)

Dict{Symbol, Float64} with 8 entries:
  :kab_ko => 148.364
  :kab_wt => 155.445
  :nb_ko  => 0.0843197
  :βr     => 0.00184341
  :r₀     => 0.806991
  :nsc₀   => 3000.0
  :nb_wt  => 0.0399102
  :pₛ     => 0.95

times = collect(0:0.1:700.0)
betamodel = ChangedBetaSelf((0.0, 700.0), times, zeros(Int64, length(times)), param_beta)

ChangedBetaSelf((0.0, 700.0), [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9  …  699.1, 699.2, 699.3, 699.4, 699.5, 699.6, 699.7, 699.8, 699.9, 700.0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0  …  0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [5 => 148.364, 4 => 155.445, 7 => 0.0843197, 3 => 0.00184341, 2 => 0.806991, 1 => 3000.0, 6 => 0.0399102, 8 => 0.95])

free_parameters(betamodel)

1-element Vector{Symbol}:
 :t₁

t₁ = 30.0
rawsim = simulate(betamodel, Dict(:t₁ => t₁))
sim = link(betamodel, rawsim)
ylabs = link_names(betamodel)
param = param_beta
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
p1 = plot(times, sim[:,1], lab="", lc=:red, lw=2, yscale=:log10, ylab=ylabs[1])
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
p2 = plot(times, sim[:,2], lab="", lc=:red, lw=2, ylab=ylabs[2])
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
p3 = plot(times, param_rt[:wt], lab="wt", lc=:black)
plot!(p3, times, rawsim[:,3], lab="del", ylab="Activation rate", lc=:red)
p4 = plot(times, hillfun(800.0).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:black)
plot!(p4, times, hillfun(t₁).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:red)
plot(p1, p2, p4, p3, xlab="Time (days)", lw=2)

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

t₁ = 400.0
rawsim = simulate(betamodel, Dict(:t₁ => t₁))
sim = link(betamodel, rawsim)
ylabs = link_names(betamodel)
param = param_beta
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
p1 = plot(times, sim[:,1], lab="", lc=:red, lw=2, yscale=:log10, ylab=ylabs[1])
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
p2 = plot(times, sim[:,2], lab="", lc=:red, lw=2, ylab=ylabs[2])
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
p3 = plot(times, param_rt[:wt], lab="wt", lc=:black)
plot!(p3, times, rawsim[:,3], lab="del", ylab="Activation rate", lc=:red)
p4 = plot(times, hillfun(800.0).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:black)
plot!(p4, times, hillfun(t₁).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:red)
plot(p1, p2, p4, p3, xlab="Time (days)", lw=2)
#savefig("example_3.svg")
savefig("betamodel_interaction.svg")

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

t₁ = 600.0
rawsim = simulate(betamodel, Dict(:t₁ => t₁))
sim = link(betamodel, rawsim)
ylabs = link_names(betamodel)
param = param_beta
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
p1 = plot(times, sim[:,1], lab="", lc=:red, lw=2, yscale=:log10, ylab=ylabs[1])
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
p2 = plot(times, sim[:,2], lab="", lc=:red, lw=2, ylab=ylabs[2])
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
p3 = plot(times, param_rt[:wt], lab="wt", lc=:black)
plot!(p3, times, rawsim[:,3], lab="del", ylab="Activation rate", lc=:red)
p4 = plot(times, hillfun(800.0).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:black)
plot!(p4, times, hillfun(t₁).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:red)
plot(p1, p2, p4, p3, xlab="Time (days)", lw=2)

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

#pgfplots()
pyplot()
itimes = collect(0:1:660)
rawsim = simulate(betamodel, Dict(:t₁ => 800.0)) # No intervention model
nonsim = link(betamodel, rawsim)
ylabs = link_names(betamodel)
param = param_beta
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
# Plot base simulation (wt)
p1 = plot(yscale=:log10, ylab=ylabs[1], xlab="Age (days)")
p2 = plot(ylab="Relative NSCs", xlab="Age (days)")
p3 = plot(xlab="Age (days)", ylab="Activation rate")
p4 = plot(xlab="Age (days)",  ylab="Self-renewal probability", legend=:bottomright)
# ???
neurogen_wt = 2 .* 0.95 .* (1 .- hillfun(800.0).(sim[:,1], times)) .* rawsim[:,2]
tot = 2 .* 0.95 .* (1 .- hillfun(800.0).(sim[:,1], times)) .* rawsim[:,2]
tot_wt = sim[:,1]
tot_rel = DataFrame(:t => times, :wt => tot_wt ./ tot_wt)
neurogen_rel = DataFrame(:t => times, :wt => neurogen_wt ./ neurogen_wt)
tot_abs = DataFrame(:t => times, :wt => tot_wt)
neurogen_abs = DataFrame(:t => times, :wt => neurogen_wt)
p5 = plot(yscale=:log10, ylab="Rate of progenitor production", xlab="Age (days)")
p6 = plot(ylab="Relative rate", xlab="Age (days)")
# Plot interventions
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
# Compute some things
step = times[2] - times[1]
non_neurogen = step .* neurogen_wt
# Preallocate vectors
loss = zeros(Float64, length(itimes))
total_neurogen = zeros(Float64, length(itimes))
for (i, itime) in enumerate(itimes)
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(itime)))
    changed = times .> itime
    sim = link(betamodel, rawsim)
    loss[i] = nonsim[findfirst(times .== 660.0), 1] .- sim[findfirst(times .== 660.0),1]
    neurogen = 2 .* 0.95 .* (1 .- hillfun(itime).(sim[:,1], times)) .* rawsim[:,2]
    total_neurogen[i] = sum(neurogen .* step) ./ sum(non_neurogen)
    plot!(p1, times[changed], sim[changed,1], lab="", lc=colours[i], la=0.9)
    plot!(p2, times, sim[:,1] ./ tot_wt, lab="", lc=colours[i], la=0.9)
    plot!(p3, times[changed], rawsim[changed,3], lab="", lc=colours[i], la=0.9)
    plot!(p4, times[changed], hillfun(itime).(sim[changed,1], times[changed]), lab="", lc=colours[i], la=0.9)
    plot!(p5, times[changed], neurogen[changed], lab="", lc=colours[i], la=0.9)
    plot!(p6, times[changed], neurogen[changed] ./ neurogen_wt[changed], lab="", lc=colours[i], la=0.9)
    tot_rel[!,"d$(Int64(round(itime)))"] = sim[:,1] ./ tot_wt
    neurogen_rel[!,"d$(Int64(round(itime)))"] = neurogen ./ neurogen_wt
    tot_abs[!,"d$(Int64(round(itime)))"] = sim[:,1]
    neurogen_abs[!,"d$(Int64(round(itime)))"] = neurogen
end
# Plot previous simulation
plot!(p1, times, nonsim[:,1], lab="", lc=:black, lw=2)
plot!(p2, times, nonsim[:,1] ./ nonsim[:,1], lab="", lc=:black, lw=2)
plot!(p3, times, rawsim[:,3], lab="", lc=:black, lw=2)
plot!(p4, times, hillfun(800.0).(nonsim[:,1], times), lab="", lc=:black, lw=2)
plot!(p5, times, neurogen_wt, lc=:black, lw=2, lab="")
plot!(p6, times, neurogen_wt ./ neurogen_wt, lc=:black, lw=2, lab="")
# Plot data
#@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
#@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
# Plot metrics
p7 = plot(itimes, loss, lab="", ylab="Stem cells lost", xlab="Intervention time (days)", lc=:black)
steps = round(minimum(total_neurogen); sigdigits=1):0.05:round(maximum(total_neurogen); sigdigits=1)
p8 = plot(itimes, total_neurogen, lc=:black, lab="", ylab="Total progenitor gain", xlab = "Intervention time (days)",
          ytick=(collect(steps), ["$(Int64(100*x))\\%" for x in steps]))
#p7 = plot(itimes, total_neurogen)
l = @layout [
    p1 p2;
    p3 p4;
    p5 p6;
    p7 p8;
]
plot(p1, p2, p4, p3, p5, p6, p7, p8, lw=1, layout=l, size=(600, 800))
#savefig("freeze_activation_individual_self.tex")
savefig("freeze_activation_individual_self.svg")
tot_rel |> CSV.write("sims/freeze_activation_individual_self_relative_nsc.csv")
neurogen_rel |> CSV.write("sims/freeze_activation_individual_self_relative_neurogen.csv")
tot_abs |> CSV.write("sims/freeze_activation_individual_self_absolute_nsc.csv")
neurogen_abs |> CSV.write("sims/freeze_activation_individual_self_absolute_neurogen.csv")

┌ Warning: No strict ticks found
└ @ PlotUtils /home/jooa/.julia/packages/PlotUtils/VgXdq/src/ticks.jl:294

"sims/freeze_activation_individual_self_absolute_neurogen.csv"

itimes = metrics_loss.itime
int_neurogenesis_1 = zeros(Float64, length(itimes))
#int_neurogenesis_2 = zeros(Float64, length(itimes))
int_neurogenesis_wt = zeros(Float64, length(itimes))
step = times[2] - times[1]
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
for (i, t₁) in enumerate(itimes)
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(t₁)))
    int_neurogenesis_1[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- hillfun(t₁).(rawsim[:,1], times)))
    #rawsim = simulate(betamodel, Dict(:t₁ => Float64(t₁)))
    #int_neurogenesis_2[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- param_bt[:wt].(times)))
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(700)))
    int_neurogenesis_wt[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- hillfun(t₁).(rawsim[:,1], times)))
end
yticks = 50:5:115
plot(xlab="Intervention Time (days)", ylab="increase in TAPs produced over wt", title="Life-long (700 days) neurogenesis", yticks=(collect(yticks), ["$x\\%" for x in yticks]))
metrics_neurogen[!,"beta_fullself"] = int_neurogenesis_1 ./ int_neurogenesis_wt
plot!(itimes, int_neurogenesis_1 ./ int_neurogenesis_wt .* 100, lab="'beta' model")
#plot!(itimes, int_neurogenesis_2 ./ int_neurogenesis_wt .* 100, lab="'value' model")
#plot!(itimes, int_neurogenesis_wt .* step, lab="wt", lc=:gray)

itimes = collect(0:10:699)#[100, 200]#, 60, 100, 300, 600]
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
rawsim = simulate(betamodel, Dict(:t₁ => 700.0))
sim = link(betamodel, rawsim)
ylabs = link_names(betamodel)
p1 = plot(times, sim[:,1], lab="", lc=:gray, lw=2, yscale=:log10, ylab=ylabs[1])
p2 = plot(times, sim[:,2], lab="", lc=:gray, lw=2, ylab=ylabs[2])
p3 = plot(times, param_rt[:wt], lab="", ylab="Activation rate", lc=:gray, lw=2)
#p4 = plot(times, param_bt[:wt], lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:black, lw=2)
p4 = plot(times, hillfun(800).(sim[:,1], times), lab="", 
          ylab="Self-renewal probability", legend=:bottomright, lc=:black, lw=2)
for (i, t₁) in enumerate(itimes)
    step = times[2] - times[1]
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(t₁)))
    sim = link(betamodel, rawsim)
    has_changed = times .> t₁
    has_changed[max(1, findfirst(has_changed)-2):end] .= true
    plot!(p1, times[has_changed], sim[has_changed,1], lab="", lw=1, lc=colours[i]) #t₁=$t₁
    plot!(p2, times[has_changed], sim[has_changed,2], lab="", lw=1, lc=colours[i])
    plot!(p3, times[has_changed], rawsim[has_changed,3], lab="", lw=1, lc=colours[i])
    plot!(p4, times[has_changed], hillfun(t₁).(sim[:,1], times)[has_changed], lab="", lw=1, lc=colours[i])
end
#plot!(p1, times, sim[:,1], lab="", lc=:red, lw=2)
@df @where(data, :genotype .== "wt", :name .== "total")  scatter!(p1, :age, :value, mc=:black, lab="")
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
plot(p1, p2, p4, p3, xlab="Time (days)", size=(600, 600))

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

#itimes = collect(0:1:660)#[100, 200]#, 60, 100, 300, 600]
itimes = metrics_loss.itime
rawsim = simulate(betamodel, Dict(:t₁ => 700.0))
regsim = link(betamodel, rawsim)
nsc660 = regsim[times .== 660.0, 1][1]
nscs = map(itimes) do t₁
    rawsim = simulate(betamodel, Dict(:t₁ => t₁))
    regsim = link(betamodel, rawsim)
    nsc660 .- regsim[times .== 660.0, 1][1]
end
metrics_loss[!,"beta_fullself"] = nscs
plot(itimes, nscs, lc=:black, lab="", xlab="Intervention Time (days)", ylab="Stem cells lost", title="Loss of stem cells from intervention")

metrics_loss

661 rows × 3 columns

	itime	type	beta_fullself
	Int64	String	Float64
1	0	aged loss	15.3519
2	1	aged loss	15.3514
3	2	aged loss	15.3514
4	3	aged loss	15.3517
5	4	aged loss	15.3517
6	5	aged loss	15.3518
7	6	aged loss	15.3519
8	7	aged loss	15.3514
9	8	aged loss	15.3515
10	9	aged loss	15.3518
11	10	aged loss	15.351
12	11	aged loss	15.3511
13	12	aged loss	15.351
14	13	aged loss	15.3514
15	14	aged loss	15.3514
16	15	aged loss	15.3513
17	16	aged loss	15.3515
18	17	aged loss	15.3519
19	18	aged loss	15.3519
20	19	aged loss	15.3523
21	20	aged loss	15.3523
22	21	aged loss	15.3517
23	22	aged loss	15.3517
24	23	aged loss	15.3517
25	24	aged loss	15.3509
26	25	aged loss	15.3508
27	26	aged loss	15.3505
28	27	aged loss	15.3504
29	28	aged loss	15.3504
30	29	aged loss	15.3512
⋮	⋮	⋮	⋮

itimes = collect(0:10:600)#[100, 200]#, 60, 100, 300, 600]
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
rawsim = simulate(betamodel, Dict(:t₁ => 700.0))
sim = link(betamodel, rawsim)
p = plot(xlab="Time (days)", ylab="Rate of TAP production from NSCs", yscale=:log10)
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
for (i, t₁) in enumerate(itimes)
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(t₁)))
    #has_changed = times .> t₁
    plot!(p, times, rawsim[:,2] .* 2 .* 0.95 .* (1 .- hillfun(t₁).(rawsim[:,1], times)), lab="", lc=colours[i])
end
rawsim = simulate(betamodel, Dict(:t₁ => Float64(700)))
plot!(p, times, rawsim[:,2] .* 2 .* 0.95 .* (1 .- param_bt[:wt].(times)), lab="WT", lc=:gray, lw=2)
#plot!(p, sim_ifnko_both.t, sim_ifnko_both.active .* sim_ifnko_both.counts .* sim_ifnko_both.b .* 2 .* 0.95, lab="KO", lc=:red, lw=2)
p
#savefig("plots/scenario1-lifelong-neurogenesis.svg")

itimes = collect(0:10:600)#[100, 200]#, 60, 100, 300, 600]
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
rawsim = simulate(betamodel, Dict(:t₁ => 700.0))
sim = link(betamodel, rawsim)
p = plot(xlab="Time (days)", ylab="Relative progenitor production", title="Progenitor production compared to wildtype")#, yscale=:log10)
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
rawsim = simulate(betamodel, Dict(:t₁ => Float64(700)))
neurogen = rawsim[:,2] .* 2 .* 0.95 .* (1 .- param_bt[:wt].(times))
for (i, t₁) in enumerate(itimes)
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(t₁)))
    #has_changed = times .> t₁
    plot!(p, times, rawsim[:,2] .* 2 .* 0.95 .* (1 .- hillfun(t₁).(rawsim[:,1], times)) ./ neurogen, lab="", lc=colours[i], cbar=true)
end
rawsim = simulate(betamodel, Dict(:t₁ => Float64(700)))
plot!(p, times, rawsim[:,2] .* 2 .* 0.95 .* (1 .- param_bt[:wt].(times)) ./ neurogen, lab="WT", lc=:black, lw=2,
      yticks=(0.5:0.25:2.5, ["$(Int64(round(x*100)))\\%" for x in 0.5:0.25:2.5]),
      xticks=(0:100:700)
)
#plot!(p, sim_ifnko_both.t, sim_ifnko_both.active .* sim_ifnko_both.counts .* sim_ifnko_both.b .* 2 .* 0.95, lab="KO", lc=:red, lw=2)
p
#savefig("plots/scenario1-lifelong-neurogenesis.svg")

itimes = collect(0:10:699)#[100, 200]#, 60, 100, 300, 600]
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
rawsim = simulate(betamodel, Dict(:t₁ => 700.0))
sim = link(betamodel, rawsim)
tot = sim[:,1]
ylabs = link_names(betamodel)
p1 = plot(lab="", lc=:gray, lw=2, ylab=ylabs[1], xlab="Time (days)", title="Total stem cells compared to wildtype")
for (i, t₁) in enumerate(itimes)
    step = times[2] - times[1]
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(t₁)))
    sim = link(betamodel, rawsim)
    plot!(p1, times, sim[:,1] ./ tot, lab="", lw=1, lc=colours[i]) #t₁=$t₁
end
plot!(p1, times, tot ./ tot, lab="", lc=:black, lw=2,
    yticks=(0.5:0.1:1.0, ["$(Int64(x*100))\\%" for x in 0.5:0.1:1.0])
)
#@df @where(data, :genotype .== "wt", :name .== "total")  scatter!(p1, :age, :value, mc=:black, lab="")
p1

itimes = metrics_loss.itime
int_neurogenesis_1 = zeros(Float64, length(itimes))
#int_neurogenesis_2 = zeros(Float64, length(itimes))
int_neurogenesis_wt = zeros(Float64, length(itimes))
step = times[2] - times[1]
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
for (i, t₁) in enumerate(itimes)
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(t₁)))
    int_neurogenesis_1[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- hillfun(t₁).(rawsim[:,1], times)))
    #rawsim = simulate(betamodel, Dict(:t₁ => Float64(t₁)))
    #int_neurogenesis_2[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- param_bt[:wt].(times)))
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(700)))
    int_neurogenesis_wt[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- hillfun(t₁).(rawsim[:,1], times)))
end
yticks = 50:5:115
plot(xlab="Intervention Time (days)", ylab="increase in TAPs produced over wt", title="Life-long (700 days) neurogenesis", yticks=(collect(yticks), ["$x\\%" for x in yticks]))
metrics_neurogen[!,"beta_fullself"] = int_neurogenesis_1 ./ int_neurogenesis_wt
plot!(itimes, int_neurogenesis_1 ./ int_neurogenesis_wt .* 100, lab="'beta' model")
#plot!(itimes, int_neurogenesis_2 ./ int_neurogenesis_wt .* 100, lab="'value' model")
#plot!(itimes, int_neurogenesis_wt .* step, lab="wt", lc=:gray)

Self-renewal is partially genotype dependent

param_beta = Dict( # from cofit
    #:nsc₀ => 1986.0,
    :nsc₀ => 3000.0,
    #:r₀ => 0.8426,
    :r₀ => 0.839538,
    #:βr => 0.001891,
    :βr => 0.00188605,
    :kab_wt => 153.81,
    :kab_ko => 153.81,
    :nb_wt => 0.039907,
    :nb_ko => 0.0841552,
    :pₛ => 0.95,
)

Dict{Symbol, Float64} with 8 entries:
  :kab_ko => 153.81
  :kab_wt => 153.81
  :nb_ko  => 0.0841552
  :βr     => 0.00188605
  :r₀     => 0.839538
  :nsc₀   => 3000.0
  :nb_wt  => 0.039907
  :pₛ     => 0.95

times = collect(0:1:700.0)
betamodel = ChangedBetaSelf((0.0, 700.0), times, zeros(Int64, length(times)), param_beta)

ChangedBetaSelf((0.0, 700.0), [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0  …  691.0, 692.0, 693.0, 694.0, 695.0, 696.0, 697.0, 698.0, 699.0, 700.0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0  …  0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [5 => 153.81, 4 => 153.81, 7 => 0.0841552, 3 => 0.00188605, 2 => 0.839538, 1 => 3000.0, 6 => 0.039907, 8 => 0.95])

free_parameters(betamodel)

1-element Vector{Symbol}:
 :t₁

t₁ = 30.0
rawsim = simulate(betamodel, Dict(:t₁ => t₁))
sim = link(betamodel, rawsim)
ylabs = link_names(betamodel)
param = param_beta
hillfun(t₁) = (x, t) -> 1 - hill(param[:kab_wt], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
p1 = plot(times, sim[:,1], lab="", lc=:red, lw=2, yscale=:log10, ylab=ylabs[1])
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
p2 = plot(times, sim[:,2], lab="", lc=:red, lw=2, ylab=ylabs[2])
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
p3 = plot(times, param_rt[:wt], lab="wt", lc=:black)
plot!(p3, times, rawsim[:,3], lab="del", ylab="Activation rate", lc=:red)
p4 = plot(times, hillfun(800.0).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:black)
plot!(p4, times, hillfun(t₁).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:red)
plot(p1, p2, p4, p3, xlab="Time (days)", lw=2)

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

t₁ = 400.0
rawsim = simulate(betamodel, Dict(:t₁ => t₁))
sim = link(betamodel, rawsim)
ylabs = link_names(betamodel)
param = param_beta
hillfun(t₁) = (x, t) -> 1 - hill(param[:kab_wt], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
p1 = plot(times, sim[:,1], lab="", lc=:red, lw=2, yscale=:log10, ylab=ylabs[1])
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
p2 = plot(times, sim[:,2], lab="", lc=:red, lw=2, ylab=ylabs[2])
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
p3 = plot(times, param_rt[:wt], lab="wt", lc=:black)
plot!(p3, times, rawsim[:,3], lab="del", ylab="Activation rate", lc=:red)
p4 = plot(times, hillfun(800.0).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:black)
plot!(p4, times, hillfun(t₁).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:red)
plot(p1, p2, p4, p3, xlab="Time (days)", lw=2)
savefig("betamodel_static.svg")

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

t₁ = 600.0
rawsim = simulate(betamodel, Dict(:t₁ => t₁))
sim = link(betamodel, rawsim)
ylabs = link_names(betamodel)
param = param_beta
hillfun(t₁) = (x, t) -> 1 - hill(param[:kab_wt], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
p1 = plot(times, sim[:,1], lab="", lc=:red, lw=2, yscale=:log10, ylab=ylabs[1])
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
p2 = plot(times, sim[:,2], lab="", lc=:red, lw=2, ylab=ylabs[2])
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
p3 = plot(times, param_rt[:wt], lab="wt", lc=:black)
plot!(p3, times, rawsim[:,3], lab="del", ylab="Activation rate", lc=:red)
p4 = plot(times, hillfun(800.0).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:black)
plot!(p4, times, hillfun(t₁).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:red)
plot(p1, p2, p4, p3, xlab="Time (days)", lw=2)

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

itimes = collect(0:10:699)#[100, 200]#, 60, 100, 300, 600]
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
rawsim = simulate(betamodel, Dict(:t₁ => 700.0))
sim = link(betamodel, rawsim)
ylabs = link_names(betamodel)
p1 = plot(times, sim[:,1], lab="", lc=:gray, lw=2, yscale=:log10, ylab=ylabs[1])
p2 = plot(times, sim[:,2], lab="", lc=:gray, lw=2, ylab=ylabs[2])
p3 = plot(times, param_rt[:wt], lab="", ylab="Activation rate", lc=:gray, lw=2)
#p4 = plot(times, param_bt[:wt], lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:black, lw=2)
p4 = plot(times, hillfun(800).(sim[:,1], times), lab="", 
          ylab="Self-renewal probability", legend=:bottomright, lc=:black, lw=2)
for (i, t₁) in enumerate(itimes)
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(t₁)))
    sim = link(betamodel, rawsim)
    has_changed = times .> t₁
    has_changed[max(1, findfirst(has_changed)-2):end] .= true
    plot!(p1, times[has_changed], sim[has_changed,1], lab="", lw=1, lc=colours[i]) #t₁=$t₁
    plot!(p2, times[has_changed], sim[has_changed,2], lab="", lw=1, lc=colours[i])
    plot!(p3, times[has_changed], rawsim[has_changed,3], lab="", lw=1, lc=colours[i])
    plot!(p4, times[has_changed], hillfun(t₁).(sim[:,1], times)[has_changed], lab="", lw=1, lc=colours[i])
end
#plot!(p1, times, sim[:,1], lab="", lc=:red, lw=2)
@df @where(data, :genotype .== "wt", :name .== "total")  scatter!(p1, :age, :value, mc=:black, lab="")
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
plot(p1, p2, p4, p3, xlab="Time (days)", size=(600, 600))

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

#pgfplots()
pyplot()
itimes = collect(0:1:660)
rawsim = simulate(betamodel, Dict(:t₁ => 800.0)) # No intervention model
nonsim = link(betamodel, rawsim)
sim = nonsim
ylabs = link_names(betamodel)
param = param_beta
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
# Plot base simulation (wt)
p1 = plot(yscale=:log10, ylab=ylabs[1], xlab="Age (days)")
p2 = plot(ylab="Relative NSCs", xlab="Age (days)")
p3 = plot(xlab="Age (days)", ylab="Activation rate")
p4 = plot(xlab="Age (days)",  ylab="Self-renewal probability", legend=:bottomright)
# ???
neurogen_wt = 2 .* 0.95 .* (1 .- hillfun(800.0).(sim[:,1], times)) .* rawsim[:,2]
tot = 2 .* 0.95 .* (1 .- hillfun(800.0).(sim[:,1], times)) .* rawsim[:,2]
tot_wt = sim[:,1]
tot_rel = DataFrame(:t => times, :wt => tot_wt ./ tot_wt)
neurogen_rel = DataFrame(:t => times, :wt => neurogen_wt ./ neurogen_wt)
tot_abs = DataFrame(:t => times, :wt => tot_wt)
neurogen_abs = DataFrame(:t => times, :wt => neurogen_wt)
p5 = plot(yscale=:log10, ylab="Rate of progenitor production", xlab="Age (days)")
p6 = plot(ylab="Relative rate", xlab="Age (days)")
# Plot interventions
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
# Compute some things
step = times[2] - times[1]
non_neurogen = step .* neurogen_wt
# Preallocate vectors
loss = zeros(Float64, length(itimes))
total_neurogen = zeros(Float64, length(itimes))
for (i, itime) in enumerate(itimes)
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(itime)))
    changed = times .> itime
    sim = link(betamodel, rawsim)
    loss[i] = nonsim[findfirst(times .== 660.0), 1] .- sim[findfirst(times .== 660.0),1]
    neurogen = 2 .* 0.95 .* (1 .- hillfun(itime).(sim[:,1], times)) .* rawsim[:,2]
    total_neurogen[i] = sum(neurogen .* step) ./ sum(non_neurogen)
    plot!(p1, times[changed], sim[changed,1], lab="", lc=colours[i], la=0.9)
    plot!(p2, times, sim[:,1] ./ tot_wt, lab="", lc=colours[i], la=0.9)
    plot!(p3, times[changed], rawsim[changed,3], lab="", lc=colours[i], la=0.9)
    plot!(p4, times[changed], hillfun(itime).(sim[changed,1], times[changed]), lab="", lc=colours[i], la=0.9)
    plot!(p5, times[changed], neurogen[changed], lab="", lc=colours[i], la=0.9)
    plot!(p6, times[changed], neurogen[changed] ./ neurogen_wt[changed], lab="", lc=colours[i], la=0.9)
    tot_rel[!,"d$(Int64(round(itime)))"] = sim[:,1] ./ tot_wt
    neurogen_rel[!,"d$(Int64(round(itime)))"] = neurogen ./ neurogen_wt
    tot_abs[!,"d$(Int64(round(itime)))"] = sim[:,1]
    neurogen_abs[!,"d$(Int64(round(itime)))"] = neurogen
end
# Plot previous simulation
plot!(p1, times, nonsim[:,1], lab="", lc=:black, lw=2)
plot!(p2, times, nonsim[:,1] ./ nonsim[:,1], lab="", lc=:black, lw=2)
plot!(p3, times, rawsim[:,3], lab="", lc=:black, lw=2)
plot!(p4, times, hillfun(800.0).(nonsim[:,1], times), lab="", lc=:black, lw=2)
plot!(p5, times, neurogen_wt, lc=:black, lw=2, lab="")
plot!(p6, times, neurogen_wt ./ neurogen_wt, lc=:black, lw=2, lab="")
# Plot data
#@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
#@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
# Plot metrics
p7 = plot(itimes, loss, lab="", ylab="Stem cells lost", xlab="Intervention time (days)", lc=:black)
steps = round(minimum(total_neurogen); sigdigits=1):0.05:round(maximum(total_neurogen); sigdigits=1)
p8 = plot(itimes, total_neurogen, lc=:black, lab="", ylab="Total progenitor gain", xlab = "Intervention time (days)",
          ytick=(collect(steps), ["$(Int64(100*x))\\%" for x in steps]))
#p7 = plot(itimes, total_neurogen)
l = @layout [
    p1 p2;
    p3 p4;
    p5 p6;
    p7 p8;
]
plot(p1, p2, p4, p3, p5, p6, p7, p8, lw=1, layout=l, size=(600, 800))
#savefig("freeze_activation_semishared_self.tex")
savefig("freeze_activation_semishared_self.svg")
tot_rel |> CSV.write("sims/freeze_activation_semishared_self_relative_nsc.csv")
neurogen_rel |> CSV.write("sims/freeze_activation_semishared_self_relative_neurogen.csv")
tot_abs |> CSV.write("sims/freeze_activation_semishared_self_absolute_nsc.csv")
neurogen_abs |> CSV.write("sims/freeze_activation_semishared_self_absolute_neurogen.csv")

┌ Warning: No strict ticks found
└ @ PlotUtils /home/jooa/.julia/packages/PlotUtils/VgXdq/src/ticks.jl:294

"sims/freeze_activation_semishared_self_absolute_neurogen.csv"

#pgfplots()
pyplot()
rawsim = simulate(betamodel, Dict(:t₁ => 800.0)) # No intervention model
itimes = collect(0:10:660)
nonsim = link(betamodel, rawsim)
ylabs = link_names(betamodel)
param = param_beta
hillfun(t₁) = (x, t) -> 1 - hill(param[:kab_wt], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
# Plot base simulation (wt)
p1 = plot(yscale=:log10, ylab=ylabs[1], xlab="Age (days)")
p2 = plot(ylab=ylabs[2], xlab="Age (days)", ylim=(0.18, 0.65))
p3 = plot(xlab="Age (days)", ylab="Activation rate")
p4 = plot(xlab="Age (days)",  ylab="Self-renewal probability", legend=:bottomright)
neurogen_wt = 2 .* 0.95 .* (1 .- hillfun(800.0).(sim[:,1], times)) .* rawsim[:,2]
tot = 2 .* 0.95 .* (1 .- hillfun(800.0).(sim[:,1], times)) .* rawsim[:,2]
tot_wt = sim[:,1]
tot_rel = DataFrame(:t => times, :wt => tot_wt ./ tot_wt)
neurogen_rel = DataFrame(:t => times, :wt => neurogen_wt ./ neurogen_wt)
tot_abs = DataFrame(:t => times, :wt => tot_wt)
neurogen_abs = DataFrame(:t => times, :wt => neurogen_wt)
# ???
neurogen = 2 .* 0.95 .* (1 .- hillfun(800.0).(sim[:,1], times)) .* rawsim[:,2]
p5 = plot(yscale=:log10, ylab="Rate of progenitor production", xlab="Age (days)")
# Plot interventions
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
# Compute some things
step = times[2] - times[1]
non_neurogen = step .* neurogen
# Preallocate vectors
loss = zeros(Float64, length(itimes))
total_neurogen = zeros(Float64, length(itimes))
for (i, itime) in enumerate(itimes)
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(itime)))
    changed = times .> itime
    sim = link(betamodel, rawsim)
    loss[i] = nonsim[findfirst(times .== 660.0), 1] .- sim[findfirst(times .== 660.0),1]
    neurogen = 2 .* 0.95 .* (1 .- hillfun(itime).(sim[:,1], times)) .* rawsim[:,2]
    total_neurogen[i] = sum(neurogen .* step) ./ sum(non_neurogen)
    plot!(p1, times[changed], sim[changed,1], lab="", lc=colours[i], la=0.9)
    plot!(p2, times, sim[:,1] ./ tot_wt, lab="", lc=colours[i], la=0.9)
    plot!(p2, times[changed], sim[changed,2], lab="", lc=colours[i], la=0.9)
    plot!(p3, times[changed], rawsim[changed,3], lab="", lc=colours[i], la=0.9)
    plot!(p4, times[changed], hillfun(itime).(sim[changed,1], times[changed]), lab="", lc=colours[i], la=0.9)
    plot!(p5, times[changed], neurogen[changed], lab="", lc=colours[i], la=0.9)
end
# Plot previous sims
plot!(p1, times, nonsim[:,1], lab="", lc=:black, lw=2)
plot!(p2, times, nonsim[:,2], lab="", lc=:black, lw=2)
plot!(p3, times, rawsim[:,3], lab="", lc=:black, lw=2)
plot!(p4, times, hillfun(800.0).(nonsim[:,1], times), lab="", lc=:black, lw=2)
plot!(p5, times, neurogen, lc=:black, lw=2, lab="")
# Plot data
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
# Plot metrics
p6 = plot(itimes, loss, lab="", ylab="Stem cells lost", xlab="Intervention time (days)", lc=:black)
neurogen_range = maximum(total_neurogen) - minimum(total_neurogen)
steps = 0.95:0.05:1.2
p7 = plot(itimes, total_neurogen, lc=:black, lab="", ylab="Total progenitor gain", xlab = "Intervention time (days)",
          ytick=(collect(steps), ["$(Int64(round(100*x)))\\%" for x in steps]))
#p7 = plot(itimes, total_neurogen)
l = @layout [
    p1 p2;
    p3 p4;
    p5{.5h} [p6;
             p7]
]
plot(p1, p2, p4, p3, p5, p6, p7, lw=1, layout=l, size=(600, 800))
#savefig("freeze_activation_individual_self.tex")
savefig("freeze_activation_semishared_self.svg")

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

itimes = collect(0:1:660)#[100, 200]#, 60, 100, 300, 600]
itimes = metrics_loss.itime
rawsim = simulate(betamodel, Dict(:t₁ => 700.0))
regsim = link(betamodel, rawsim)
nsc660 = regsim[times .== 660.0, 1][1]
nscs = map(itimes) do t₁
    rawsim = simulate(betamodel, Dict(:t₁ => t₁))
    regsim = link(betamodel, rawsim)
    nsc660 .- regsim[times .== 660.0, 1][1]
end
metrics_loss[!,"beta_partself"] = nscs
plot(itimes, nscs, lc=:black, lab="", xlab="Intervention Time (days)", ylab="Stem cells lost", title="Loss of stem cells from intervention")

itimes = collect(0:10:600)#[100, 200]#, 60, 100, 300, 600]
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
rawsim = simulate(betamodel, Dict(:t₁ => 700.0))
sim = link(betamodel, rawsim)
p = plot(xlab="Time (days)", ylab="Rate of TAP production from NSCs", yscale=:log10)
hillfun(t₁) = (x, t) -> 1 - hill(param[:kab_wt], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
for (i, t₁) in enumerate(itimes)
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(t₁)))
    #has_changed = times .> t₁
    plot!(p, times, rawsim[:,2] .* 2 .* 0.95 .* (1 .- hillfun(t₁).(rawsim[:,1], times)), lab="", lc=colours[i])
end
rawsim = simulate(betamodel, Dict(:t₁ => Float64(700)))
plot!(p, times, rawsim[:,2] .* 2 .* 0.95 .* (1 .- param_bt[:wt].(times)), lab="WT", lc=:gray, lw=2)
#plot!(p, sim_ifnko_both.t, sim_ifnko_both.active .* sim_ifnko_both.counts .* sim_ifnko_both.b .* 2 .* 0.95, lab="KO", lc=:red, lw=2)
p
#savefig("plots/scenario1-lifelong-neurogenesis.svg")

itimes = metrics_loss.itime
int_neurogenesis_1 = zeros(Float64, length(itimes))
#int_neurogenesis_2 = zeros(Float64, length(itimes))
int_neurogenesis_wt = zeros(Float64, length(itimes))
step = times[2] - times[1]
hillfun(t₁) = (x, t) -> 1 - hill(param[:kab_wt], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
for (i, t₁) in enumerate(itimes)
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(t₁)))
    int_neurogenesis_1[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- hillfun(t₁).(rawsim[:,1], times)))
    #rawsim = simulate(betamodel, Dict(:t₁ => Float64(t₁)))
    #int_neurogenesis_2[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- param_bt[:wt].(times)))
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(700)))
    int_neurogenesis_wt[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- hillfun(t₁).(rawsim[:,1], times)))
end
yticks = 50:5:115
plot(xlab="Intervention Time (days)", ylab="increase in TAPs produced over wt", title="Life-long (700 days) neurogenesis", yticks=(collect(yticks), ["$x\\%" for x in yticks]))
metrics_neurogen[!,"beta_partself"] = int_neurogenesis_1 ./ int_neurogenesis_wt
plot!(itimes, int_neurogenesis_1 ./ int_neurogenesis_wt .* 100, lab="'beta' model")
#plot!(itimes, int_neurogenesis_2 ./ int_neurogenesis_wt .* 100, lab="'value' model")
#plot!(itimes, int_neurogenesis_wt .* step, lab="wt", lc=:gray)

Self-renewal is genotype invariant

param_beta = Dict( # from cofit
    #:nsc₀ => 1986.0,
    :nsc₀ => 5000.0,
    #:r₀ => 0.8426,
    :r₀ => 0.672387,
    #:βr => 0.001891,
    :βr => 0.0015538,
    :kab_wt => 153.724,
    :kab_ko => 153.724,
    :nb_wt => 0.0584954,
    :nb_ko => 0.0584954,
    #:nb_ko => 0.0765863,
    :pₛ => 0.95,
)

Dict{Symbol, Float64} with 8 entries:
  :kab_ko => 153.724
  :kab_wt => 153.724
  :nb_ko  => 0.0584954
  :βr     => 0.0015538
  :r₀     => 0.672387
  :nsc₀   => 5000.0
  :nb_wt  => 0.0584954
  :pₛ     => 0.95

times = collect(0:1:700.0)
betamodel = ChangedBetaSelf((0.0, 700.0), times, zeros(Int64, length(times)), param_beta)

ChangedBetaSelf((0.0, 700.0), [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0  …  691.0, 692.0, 693.0, 694.0, 695.0, 696.0, 697.0, 698.0, 699.0, 700.0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0  …  0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [5 => 153.724, 4 => 153.724, 7 => 0.0584954, 3 => 0.0015538, 2 => 0.672387, 1 => 5000.0, 6 => 0.0584954, 8 => 0.95])

free_parameters(betamodel)

1-element Vector{Symbol}:
 :t₁

t₁ = 30.0
rawsim = simulate(betamodel, Dict(:t₁ => t₁))
sim = link(betamodel, rawsim)
ylabs = link_names(betamodel)
param = param_beta
hillfun(t₁) = (x, t) -> 1 - hill(param[:kab_wt], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
p1 = plot(times, sim[:,1], lab="", lc=:red, lw=2, yscale=:log10, ylab=ylabs[1])
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
p2 = plot(times, sim[:,2], lab="", lc=:red, lw=2, ylab=ylabs[2])
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
p3 = plot(times, param_rt[:wt], lab="wt", lc=:black)
plot!(p3, times, rawsim[:,3], lab="del", ylab="Activation rate", lc=:red)
p4 = plot(times, hillfun(800.0).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:black)
plot!(p4, times, hillfun(t₁).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:red)
plot(p1, p2, p4, p3, xlab="Time (days)", lw=2)

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

t₁ = 400.0
rawsim = simulate(betamodel, Dict(:t₁ => t₁))
sim = link(betamodel, rawsim)
ylabs = link_names(betamodel)
param = param_beta
hillfun(t₁) = (x, t) -> 1 - hill(param[:kab_wt], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
p1 = plot(times, sim[:,1], lab="", lc=:red, lw=2, yscale=:log10, ylab=ylabs[1])
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
p2 = plot(times, sim[:,2], lab="", lc=:red, lw=2, ylab=ylabs[2])
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
p3 = plot(times, param_rt[:wt], lab="wt", lc=:black)
plot!(p3, times, rawsim[:,3], lab="del", ylab="Activation rate", lc=:red)
p4 = plot(times, hillfun(800.0).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:black)
plot!(p4, times, hillfun(t₁).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:red)
plot(p1, p2, p4, p3, xlab="Time (days)", lw=2)
savefig("betamodel_static.svg")

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

t₁ = 600.0
rawsim = simulate(betamodel, Dict(:t₁ => t₁))
sim = link(betamodel, rawsim)
ylabs = link_names(betamodel)
param = param_beta
hillfun(t₁) = (x, t) -> 1 - hill(param[:kab_wt], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
p1 = plot(times, sim[:,1], lab="", lc=:red, lw=2, yscale=:log10, ylab=ylabs[1])
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
p2 = plot(times, sim[:,2], lab="", lc=:red, lw=2, ylab=ylabs[2])
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
p3 = plot(times, param_rt[:wt], lab="wt", lc=:black)
plot!(p3, times, rawsim[:,3], lab="del", ylab="Activation rate", lc=:red)
p4 = plot(times, hillfun(800.0).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:black)
plot!(p4, times, hillfun(t₁).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:red)
plot(p1, p2, p4, p3, xlab="Time (days)", lw=2)

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

itimes = collect(0:10:699)#[100, 200]#, 60, 100, 300, 600]
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
rawsim = simulate(betamodel, Dict(:t₁ => 700.0))
sim = link(betamodel, rawsim)
ylabs = link_names(betamodel)
p1 = plot(times, sim[:,1], lab="", lc=:gray, lw=2, yscale=:log10, ylab=ylabs[1])
p2 = plot(times, sim[:,2], lab="", lc=:gray, lw=2, ylab=ylabs[2])
p3 = plot(times, param_rt[:wt], lab="", ylab="Activation rate", lc=:gray, lw=2)
#p4 = plot(times, param_bt[:wt], lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:black, lw=2)
p4 = plot(times, hillfun(800).(sim[:,1], times), lab="", 
          ylab="Self-renewal probability", legend=:bottomright, lc=:black, lw=2)
for (i, t₁) in enumerate(itimes)
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(t₁)))
    sim = link(betamodel, rawsim)
    has_changed = times .> t₁
    has_changed[max(1, findfirst(has_changed)-2):end] .= true
    plot!(p1, times[has_changed], sim[has_changed,1], lab="", lw=1, lc=colours[i]) #t₁=$t₁
    plot!(p2, times[has_changed], sim[has_changed,2], lab="", lw=1, lc=colours[i])
    plot!(p3, times[has_changed], rawsim[has_changed,3], lab="", lw=1, lc=colours[i])
    plot!(p4, times[has_changed], hillfun(t₁).(sim[:,1], times)[has_changed], lab="", lw=1, lc=colours[i])
end
#plot!(p1, times, sim[:,1], lab="", lc=:red, lw=2)
@df @where(data, :genotype .== "wt", :name .== "total")  scatter!(p1, :age, :value, mc=:black, lab="")
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
plot(p1, p2, p4, p3, xlab="Time (days)", size=(600, 600))

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

#pgfplots()
pyplot()
itimes = collect(0:1:660)
rawsim = simulate(betamodel, Dict(:t₁ => 800.0)) # No intervention model
nonsim = link(betamodel, rawsim)
sim = nonsim
ylabs = link_names(betamodel)
param = param_beta
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
# Plot base simulation (wt)
p1 = plot(yscale=:log10, ylab=ylabs[1], xlab="Age (days)")
p2 = plot(ylab="Relative NSCs", xlab="Age (days)")
p3 = plot(xlab="Age (days)", ylab="Activation rate")
p4 = plot(xlab="Age (days)",  ylab="Self-renewal probability", legend=:bottomright)
# ???
neurogen_wt = 2 .* 0.95 .* (1 .- hillfun(800.0).(sim[:,1], times)) .* rawsim[:,2]
tot = 2 .* 0.95 .* (1 .- hillfun(800.0).(sim[:,1], times)) .* rawsim[:,2]
tot_wt = sim[:,1]
tot_rel = DataFrame(:t => times, :wt => tot_wt ./ tot_wt)
neurogen_rel = DataFrame(:t => times, :wt => neurogen_wt ./ neurogen_wt)
tot_abs = DataFrame(:t => times, :wt => tot_wt)
neurogen_abs = DataFrame(:t => times, :wt => neurogen_wt)
p5 = plot(yscale=:log10, ylab="Rate of progenitor production", xlab="Age (days)")
p6 = plot(ylab="Relative rate", xlab="Age (days)")
# Plot interventions
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
# Compute some things
step = times[2] - times[1]
non_neurogen = step .* neurogen_wt
# Preallocate vectors
loss = zeros(Float64, length(itimes))
total_neurogen = zeros(Float64, length(itimes))
for (i, itime) in enumerate(itimes)
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(itime)))
    changed = times .> itime
    sim = link(betamodel, rawsim)
    loss[i] = nonsim[findfirst(times .== 660.0), 1] .- sim[findfirst(times .== 660.0),1]
    neurogen = 2 .* 0.95 .* (1 .- hillfun(itime).(sim[:,1], times)) .* rawsim[:,2]
    total_neurogen[i] = sum(neurogen .* step) ./ sum(non_neurogen)
    plot!(p1, times[changed], sim[changed,1], lab="", lc=colours[i], la=0.9)
    plot!(p2, times, sim[:,1] ./ tot_wt, lab="", lc=colours[i], la=0.9)
    plot!(p3, times[changed], rawsim[changed,3], lab="", lc=colours[i], la=0.9)
    plot!(p4, times[changed], hillfun(itime).(sim[changed,1], times[changed]), lab="", lc=colours[i], la=0.9)
    plot!(p5, times[changed], neurogen[changed], lab="", lc=colours[i], la=0.9)
    plot!(p6, times[changed], neurogen[changed] ./ neurogen_wt[changed], lab="", lc=colours[i], la=0.9)
    tot_rel[!,"d$(Int64(round(itime)))"] = sim[:,1] ./ tot_wt
    neurogen_rel[!,"d$(Int64(round(itime)))"] = neurogen ./ neurogen_wt
    tot_abs[!,"d$(Int64(round(itime)))"] = sim[:,1]
    neurogen_abs[!,"d$(Int64(round(itime)))"] = neurogen
end
# Plot previous simulation
plot!(p1, times, nonsim[:,1], lab="", lc=:black, lw=2)
plot!(p2, times, nonsim[:,1] ./ nonsim[:,1], lab="", lc=:black, lw=2)
plot!(p3, times, rawsim[:,3], lab="", lc=:black, lw=2)
plot!(p4, times, hillfun(800.0).(nonsim[:,1], times), lab="", lc=:black, lw=2)
plot!(p5, times, neurogen_wt, lc=:black, lw=2, lab="")
plot!(p6, times, neurogen_wt ./ neurogen_wt, lc=:black, lw=2, lab="")
# Plot data
#@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
#@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
# Plot metrics
p7 = plot(itimes, loss, lab="", ylab="Stem cells lost", xlab="Intervention time (days)", lc=:black)
steps = round(minimum(total_neurogen); sigdigits=1):0.05:round(maximum(total_neurogen); sigdigits=1)
p8 = plot(itimes, total_neurogen, lc=:black, lab="", ylab="Total progenitor gain", xlab = "Intervention time (days)",
          ytick=(collect(steps), ["$(Int64(100*x))\\%" for x in steps]))
#p7 = plot(itimes, total_neurogen)
l = @layout [
    p1 p2;
    p3 p4;
    p5 p6;
    p7 p8;
]
plot(p1, p2, p4, p3, p5, p6, p7, p8, lw=1, layout=l, size=(600, 800))
#savefig("freeze_activation_shared_self.tex")
savefig("freeze_activation_shared_self.svg")
tot_rel |> CSV.write("sims/freeze_activation_shared_self_relative_nsc.csv")
neurogen_rel |> CSV.write("sims/freeze_activation_shared_self_relative_neurogen.csv")
tot_abs |> CSV.write("sims/freeze_activation_shared_self_absolute_nsc.csv")
neurogen_abs |> CSV.write("sims/freeze_activation_shared_self_absolute_neurogen.csv")

┌ Warning: No strict ticks found
└ @ PlotUtils /home/jooa/.julia/packages/PlotUtils/VgXdq/src/ticks.jl:294

"sims/freeze_activation_shared_self_absolute_neurogen.csv"

#pgfplots()
pyplot()
rawsim = simulate(betamodel, Dict(:t₁ => 800.0)) # No intervention model
itimes = collect(0:10:660)
nonsim = link(betamodel, rawsim)
ylabs = link_names(betamodel)
param = param_beta
hillfun(t₁) = (x, t) -> 1 - hill(param[:kab_wt], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
# Plot base simulation (wt)
p1 = plot(yscale=:log10, ylab=ylabs[1], xlab="Age (days)")
#p2 = plot(ylab=ylabs[2], xlab="Age (days)")
p2 = plot(ylab=ylabs[2], xlab="Age (days)", ylim=(0.18, 0.65))
p3 = plot(xlab="Age (days)", ylab="Activation rate")
p4 = plot(xlab="Age (days)",  ylab="Self-renewal probability", legend=:bottomright)
# ???
neurogen = 2 .* 0.95 .* (1 .- hillfun(800.0).(sim[:,1], times)) .* rawsim[:,2]
p5 = plot(yscale=:log10, ylab="Rate of progenitor production", xlab="Age (days)")
# Plot interventions
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
# Compute some things
step = times[2] - times[1]
non_neurogen = step .* neurogen
# Preallocate vectors
loss = zeros(Float64, length(itimes))
total_neurogen = zeros(Float64, length(itimes))
for (i, itime) in enumerate(itimes)
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(itime)))
    changed = times .> itime
    sim = link(betamodel, rawsim)
    loss[i] = nonsim[findfirst(times .== 660.0), 1] .- sim[findfirst(times .== 660.0),1]
    neurogen = 2 .* 0.95 .* (1 .- hillfun(itime).(sim[:,1], times)) .* rawsim[:,2]
    total_neurogen[i] = sum(neurogen .* step) ./ sum(non_neurogen)
    plot!(p1, times[changed], sim[changed,1], lab="", lc=colours[i], la=0.9)
    plot!(p2, times[changed], sim[changed,2], lab="", lc=colours[i], la=0.9)
    plot!(p3, times[changed], rawsim[changed,3], lab="", lc=colours[i], la=0.9)
    plot!(p4, times[changed], hillfun(itime).(sim[changed,1], times[changed]), lab="", lc=colours[i], la=0.9)
    plot!(p5, times[changed], neurogen[changed], lab="", lc=colours[i], la=0.9)
end
# Plot previous sims
plot!(p1, times, nonsim[:,1], lab="", lc=:black, lw=2)
plot!(p2, times, nonsim[:,2], lab="", lc=:black, lw=2)
plot!(p3, times, rawsim[:,3], lab="", lc=:black, lw=2)
plot!(p4, times, hillfun(800.0).(nonsim[:,1], times), lab="", lc=:black, lw=2)
plot!(p5, times, neurogen, lc=:black, lw=2, lab="")
# Plot data
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
# Plot metrics
p6 = plot(itimes, loss, lab="", ylab="Stem cells lost", xlab="Intervention time (days)", lc=:black)
neurogen_range = maximum(total_neurogen) - minimum(total_neurogen)
steps = 0.95:0.05:1.2
p7 = plot(itimes, total_neurogen, lc=:black, lab="", ylab="Total progenitor gain", xlab = "Intervention time (days)",
          ytick=(collect(steps), ["$(Int64(round(100*x)))\\%" for x in steps]))
#p7 = plot(itimes, total_neurogen)
l = @layout [
    p1 p2;
    p3 p4;
    p5{.5h} [p6;
             p7]
]
plot(p1, p2, p4, p3, p5, p6, p7, lw=1, layout=l, size=(600, 800))
#savefig("freeze_activation_individual_self.tex")
savefig("freeze_activation_shared_self.svg")

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

itimes = collect(0:1:660)#[100, 200]#, 60, 100, 300, 600]
itimes = metrics_loss.itime
rawsim = simulate(betamodel, Dict(:t₁ => 700.0))
regsim = link(betamodel, rawsim)
nsc660 = regsim[times .== 660.0, 1][1]
nscs = map(itimes) do t₁
    rawsim = simulate(betamodel, Dict(:t₁ => t₁))
    regsim = link(betamodel, rawsim)
    nsc660 .- regsim[times .== 660.0, 1][1]
end
metrics_loss[!,"beta_noself"] = nscs
plot(itimes, nscs, lc=:black, lab="", xlab="Intervention Time (days)", ylab="Stem cells lost", title="Loss of stem cells from intervention")

itimes = collect(0:10:600)#[100, 200]#, 60, 100, 300, 600]
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
rawsim = simulate(betamodel, Dict(:t₁ => 700.0))
sim = link(betamodel, rawsim)
p = plot(xlab="Time (days)", ylab="Rate of TAP production from NSCs", yscale=:log10)
#hillfun(t₁) = (x, t) -> 1 - hill(param[:kab], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
for (i, t₁) in enumerate(itimes)
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(t₁)))
    #has_changed = times .> t₁
    plot!(p, times, rawsim[:,2] .* 2 .* 0.95 .* (1 .- hillfun(t₁).(rawsim[:,1], times)), lab="", lc=colours[i])
end
rawsim = simulate(betamodel, Dict(:t₁ => Float64(700)))
plot!(p, times, rawsim[:,2] .* 2 .* 0.95 .* (1 .- param_bt[:wt].(times)), lab="WT", lc=:gray, lw=2)
#plot!(p, sim_ifnko_both.t, sim_ifnko_both.active .* sim_ifnko_both.counts .* sim_ifnko_both.b .* 2 .* 0.95, lab="KO", lc=:red, lw=2)
p
#savefig("plots/scenario1-lifelong-neurogenesis.svg")

itimes = metrics_loss.itime
int_neurogenesis_1 = zeros(Float64, length(itimes))
#int_neurogenesis_2 = zeros(Float64, length(itimes))
int_neurogenesis_wt = zeros(Float64, length(itimes))
step = times[2] - times[1]
#hillfun(t₁) = (x, t) -> 1 - hill(param[:kab], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
for (i, t₁) in enumerate(itimes)
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(t₁)))
    int_neurogenesis_1[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- hillfun(t₁).(rawsim[:,1], times)))
    #rawsim = simulate(betamodel, Dict(:t₁ => Float64(t₁)))
    #int_neurogenesis_2[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- param_bt[:wt].(times)))
    rawsim = simulate(betamodel, Dict(:t₁ => Float64(700)))
    int_neurogenesis_wt[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- hillfun(t₁).(rawsim[:,1], times)))
end
yticks = 50:5:115
plot(xlab="Intervention Time (days)", ylab="increase in TAPs produced over wt", title="Life-long (700 days) neurogenesis", yticks=(collect(yticks), ["$x\\%" for x in yticks]))
metrics_neurogen[!,"beta_noself"] = int_neurogenesis_1 ./ int_neurogenesis_wt
plot!(itimes, int_neurogenesis_1 ./ int_neurogenesis_wt .* 100, lab="'beta' model")
#plot!(itimes, int_neurogenesis_2 ./ int_neurogenesis_wt .* 100, lab="'value' model")
#plot!(itimes, int_neurogenesis_wt .* step, lab="wt", lc=:gray)

Value/self model

In this model the effect of interferon is a changed $r$.

This just means that we need a two-stage $r(t)$ with parameters $\theta_r = (r_0, \beta_r, r_1, t_1)$ which then is $$ r(t) = \begin{cases} r_0 \cdot e^{-\beta_rt} & \text{if } t < t_1\\ r_1 & \text{else}\\ \end{cases} $$

function broken_exponential_decay(x₀, β, x₁, t₁, t)
    if t < t₁
        x₀ * exp(-β*t)
    else
        x₁
    end
end
broken_exponential_decay(x₀, β, x₁, t₁) = t -> broken_exponential_decay(x₀, β, x₁, t₁, t)

broken_exponential_decay (generic function with 2 methods)

struct ChangedValueSelf <: ODEModel
    tspan::Tuple{Float64, Float64}
    times::Vector{Float64}
    values::Vector{Int64}
    fixed::Vector{Pair{Int64, Float64}}
    ChangedValueSelf(tspan, times, values, fixed::Vector{Pair{I, N}}) where {I <: Integer, N <: Number} = new(tspan, times, values, fixed)
end
function ratefun(model::ChangedValueSelf)
    function(du, u, p, t)
        _, r₀, βr, kab_wt, kab_ko, nb_wt, nb_ko, pₛ, r₁, t₁ = p
        Q = u[1]; A = u[2]
        b = 1 - hill(t < t₁ ? kab_wt : kab_ko, 
                     t < t₁ ? nb_wt : nb_ko, Q+A)
        r = broken_exponential_decay(r₀, βr, r₁, t₁, t)
        du[1] = dQ = -r * Q + 2* b * pₛ * A
        du[2] = dA = r * Q - pₛ * A
    end
end
function initial(t::ChangedValueSelf, x::AbstractVector)
    nsc₀, r₀, _, kab_wt, _, nb_wt, _, pₛ,  _, _ = x
    b₀ = hill(kab_wt, nb_wt, 0)
    ratio = sqrt(((pₛ - r₀)/(2*r₀))^2 + (2*b₀*pₛ) / r₀)
    nsc₀ .* [1-1/(ratio+1), 1/(ratio+1)]
end
link(t::ChangedValueSelf, x::AbstractArray) = hcat(x[:,1] .+ x[:,2], x[:,2] ./ (x[:,1] .+ x[:,2]))
parameter_names(t::Type{ChangedValueSelf}) = [:nsc₀, :r₀, :βr, :kab_wt, :kab_ko, :nb_wt, :nb_ko, :pₛ, :r₁, :t₁]
bounds(t::ChangedValueSelf) = [(100.0, 10000.0), (0.0, 1.0), (0.0, 0.1), (0.0, 1000.0), (0.0, 1000.0), (0.01, 0.1), (0.01, 0.1), (0.0, 1.0), (0.0, 1.0), (0.0, 700.0)]
output_names(t::ChangedValueSelf) = ["qNSC", "aNSC"]
link_names(t::ChangedValueSelf) = ["Total NSC", "Fraction active NSC"]

link_names (generic function with 4 methods)

Self-renewal depends on genotype

#=
param_value = Dict( # from partial cofit
    :nsc₀ => 1986.0,
    :r₀ => 0.8426,
    :r₁ => 0.4928,
    :βr => 0.001891,
    :kab => 147.969,
    :nb_wt => 0.0300083,
    :nb_ko => 0.0765863,
    :pₛ => 0.95,
)
=#
param_value = Dict( # from partial cofit
    :nsc₀ => 3000.0,
    :r₀ => 0.80691,
    :r₁ => 0.493757,
    :βr => 0.00184341,
    :kab_wt => 148.364,
    :kab_ko => 155.445,
    :nb_wt => 0.0399102,
    :nb_ko => 0.0843197,
    :pₛ => 0.95,
)

Dict{Symbol, Float64} with 9 entries:
  :kab_ko => 155.445
  :kab_wt => 148.364
  :nb_ko  => 0.0843197
  :βr     => 0.00184341
  :r₀     => 0.80691
  :nsc₀   => 3000.0
  :r₁     => 0.493757
  :nb_wt  => 0.0399102
  :pₛ     => 0.95

times = collect(0:1:700.0)
valuemodel = ChangedValueSelf((0.0, 700.0), times, zeros(Int64, length(times)), param_value)

ChangedValueSelf((0.0, 700.0), [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0  …  691.0, 692.0, 693.0, 694.0, 695.0, 696.0, 697.0, 698.0, 699.0, 700.0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0  …  0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [5 => 155.445, 4 => 148.364, 7 => 0.0843197, 3 => 0.00184341, 2 => 0.80691, 1 => 3000.0, 9 => 0.493757, 6 => 0.0399102, 8 => 0.95])

free_parameters(valuemodel)

1-element Vector{Symbol}:
 :t₁

valuemodel_r = t₁ -> broken_exponential_decay(param_value[:r₀], param_value[:βr], param_value[:r₁], t₁)

#431 (generic function with 1 method)

#pgfplots()
pyplot()
itimes = collect(0:1:660)
rawsim = simulate(valuemodel, Dict(:t₁ => 800.0)) # No intervention model
nonsim = link(valuemodel, rawsim)
sim = nonsim
ylabs = link_names(valuemodel)
param = param_beta
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
# Plot base simulation (wt)
p1 = plot(yscale=:log10, ylab=ylabs[1], xlab="Age (days)")
p2 = plot(ylab="Relative NSCs", xlab="Age (days)")
p3 = plot(xlab="Age (days)", ylab="Activation rate")
p4 = plot(xlab="Age (days)",  ylab="Self-renewal probability", legend=:bottomright)
# ???
neurogen_wt = 2 .* 0.95 .* (1 .- hillfun(800.0).(sim[:,1], times)) .* rawsim[:,2]
tot = 2 .* 0.95 .* (1 .- hillfun(800.0).(sim[:,1], times)) .* rawsim[:,2]
tot_wt = sim[:,1]
tot_rel = DataFrame(:t => times, :wt => tot_wt ./ tot_wt)
neurogen_rel = DataFrame(:t => times, :wt => neurogen_wt ./ neurogen_wt)
tot_abs = DataFrame(:t => times, :wt => tot_wt)
neurogen_abs = DataFrame(:t => times, :wt => neurogen_wt)
p5 = plot(yscale=:log10, ylab="Rate of progenitor production", xlab="Age (days)")
p6 = plot(ylab="Relative rate", xlab="Age (days)")
# Plot interventions
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
# Compute some things
step = times[2] - times[1]
non_neurogen = step .* neurogen_wt
# Preallocate vectors
loss = zeros(Float64, length(itimes))
total_neurogen = zeros(Float64, length(itimes))
for (i, itime) in enumerate(itimes)
    rawsim = simulate(valuemodel, Dict(:t₁ => Float64(itime)))
    changed = times .> itime
    sim = link(valuemodel, rawsim)
    loss[i] = nonsim[findfirst(times .== 660.0), 1] .- sim[findfirst(times .== 660.0),1]
    neurogen = 2 .* 0.95 .* (1 .- hillfun(itime).(sim[:,1], times)) .* rawsim[:,2]
    total_neurogen[i] = sum(neurogen .* step) ./ sum(non_neurogen)
    plot!(p1, times[changed], sim[changed,1], lab="", lc=colours[i], la=0.9)
    plot!(p2, times, sim[:,1] ./ tot_wt, lab="", lc=colours[i], la=0.9)
    #plot!(p3, times[changed], rawsim[changed,3], lab="", lc=colours[i], la=0.9)
    plot!(p3, times[changed], valuemodel_r(itime).(times[changed]), lab="", lc=colours[i], la=0.9)
    plot!(p4, times[changed], hillfun(itime).(sim[changed,1], times[changed]), lab="", lc=colours[i], la=0.9)
    plot!(p5, times[changed], neurogen[changed], lab="", lc=colours[i], la=0.9)
    plot!(p6, times, neurogen ./ neurogen_wt, lab="", lc=colours[i], la=0.9)
    tot_rel[!,"d$(Int64(round(itime)))"] = sim[:,1] ./ tot_wt
    neurogen_rel[!,"d$(Int64(round(itime)))"] = neurogen ./ neurogen_wt
    tot_abs[!,"d$(Int64(round(itime)))"] = sim[:,1]
    neurogen_abs[!,"d$(Int64(round(itime)))"] = neurogen
end
# Plot previous simulation
plot!(p1, times, nonsim[:,1], lab="", lc=:black, lw=2)
plot!(p2, times, nonsim[:,1] ./ nonsim[:,1], lab="", lc=:black, lw=2)
plot!(p3, times, valuemodel_r(900.0), lab="", lc=:black, lw=2)
plot!(p4, times, hillfun(800.0).(nonsim[:,1], times), lab="", lc=:black, lw=2)
plot!(p5, times, neurogen_wt, lc=:black, lw=2, lab="")
plot!(p6, times, neurogen_wt ./ neurogen_wt, lc=:black, lw=2, lab="")
# Plot data
#@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
#@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
# Plot metrics
p7 = plot(itimes, loss, lab="", ylab="Stem cells lost", xlab="Intervention time (days)", lc=:black)
steps = round(minimum(total_neurogen); sigdigits=1):0.05:round(maximum(total_neurogen); sigdigits=1)
p8 = plot(itimes, total_neurogen, lc=:black, lab="", ylab="Total progenitor gain", xlab = "Intervention time (days)",
          ytick=(collect(steps), ["$(Int64(100*x))\\%" for x in steps]))
#p7 = plot(itimes, total_neurogen)
l = @layout [
    p1 p2;
    p3 p4;
    p5 p6;
    p7 p8;
]
plot(p1, p2, p4, p3, p5, p6, p7, p8, lw=1, layout=l, size=(600, 800))
#savefig("jump_activation_individual_self.tex")
savefig("jump_activation_individual_self.svg")
tot_rel |> CSV.write("sims/jump_activation_individual_self_relative_nsc.csv")
neurogen_rel |> CSV.write("sims/jump_activation_individual_self_relative_neurogen.csv")
tot_abs |> CSV.write("sims/jump_activation_individual_self_absolute_nsc.csv")
neurogen_abs |> CSV.write("sims/jump_activation_individual_self_absolute_neurogen.csv")

┌ Warning: No strict ticks found
└ @ PlotUtils /home/jooa/.julia/packages/PlotUtils/VgXdq/src/ticks.jl:294

"sims/jump_activation_individual_self_absolute_neurogen.csv"

#pgfplots()
pyplot()
rawsim = simulate(valuemodel, Dict(:t₁ => 800.0)) # No intervention model
itimes = collect(0:10:660)
nonsim = link(valuemodel, rawsim)
ylabs = link_names(valuemodel)
param = param_beta
#hillfun(t₁) = (x, t) -> 1 - hill(param[:kab], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
# Plot base simulation (wt)
p1 = plot(yscale=:log10, ylab=ylabs[1], xlab="Age (days)")
#p2 = plot(ylab=ylabs[2], xlab="Age (days)")
p2 = plot(ylab=ylabs[2], xlab="Age (days)", ylim=(0.18, 0.65))
p3 = plot(xlab="Age (days)", ylab="Activation rate")
p4 = plot(xlab="Age (days)",  ylab="Self-renewal probability", legend=:bottomright)
# ???
neurogen = 2 .* 0.95 .* (1 .- hillfun(800.0).(sim[:,1], times)) .* rawsim[:,2]
p5 = plot(yscale=:log10, ylab="Rate of progenitor production", xlab="Age (days)")
# Plot interventions
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
# Compute some things
step = times[2] - times[1]
non_neurogen = step .* neurogen
# Preallocate vectors
loss = zeros(Float64, length(itimes))
total_neurogen = zeros(Float64, length(itimes))
for (i, itime) in enumerate(itimes)
    rawsim = simulate(valuemodel, Dict(:t₁ => Float64(itime)))
    changed = times .> itime
    changed[max(1, findfirst(changed)-20):end] .= true
    sim = link(valuemodel, rawsim)
    loss[i] = nonsim[findfirst(times .== 660.0), 1] .- sim[findfirst(times .== 660.0),1]
    neurogen = 2 .* 0.95 .* (1 .- hillfun(itime).(sim[:,1], times)) .* rawsim[:,2]
    total_neurogen[i] = sum(neurogen .* step) ./ sum(non_neurogen)
    plot!(p1, times[changed], sim[changed,1], lab="", lc=colours[i], la=0.9)
    plot!(p2, times[changed], sim[changed,2], lab="", lc=colours[i], la=0.9)
#    plot!(p3, times[changed], rawsim[changed,3], lab="", lc=colours[i], la=0.9)
    plot!(p3, times[changed], valuemodel_r(itime).(times[changed]), lab="", lc=colours[i], la=0.9)
    plot!(p4, times[changed], hillfun(itime).(sim[changed,1], times[changed]), lab="", lc=colours[i], la=0.9)
    plot!(p5, times[changed], neurogen[changed], lab="", lc=colours[i], la=0.9)
end
# Plot previous sims
plot!(p1, times, nonsim[:,1], lab="", lc=:black, lw=2)
plot!(p2, times, nonsim[:,2], lab="", lc=:black, lw=2)
#p3 = plot(times, param_rt[:wt], lab="wt", lc=:black)
#plot!(p3, times, valuemodel_r(t₁).(times), lab="del", ylab="Activation rate", lc=:red)
plot!(p3, times, valuemodel_r(900.0), lab="", lc=:black, lw=2)
plot!(p4, times, hillfun(800.0).(nonsim[:,1], times), lab="", lc=:black, lw=2)
plot!(p5, times, neurogen, lc=:black, lw=2,lab="")
# Plot data
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
# Plot metrics
p6 = plot(itimes, loss, lab="", ylab="Stem cells lost", xlab="Intervention time (days)", lc=:black)
neurogen_range = maximum(total_neurogen) - minimum(total_neurogen)
steps = 0.70:0.1:1.1
p7 = plot(itimes, total_neurogen, lc=:black, lab="", ylab="Total progenitor gain", xlab = "Intervention time (days)",
          ytick=(collect(steps), ["$(Int64(round(100*x)))\\%" for x in steps]))
#p7 = plot(itimes, total_neurogen)
l = @layout [
    p1 p2;
    p3 p4;
    p5{.5h} [p6;
             p7]
]
plot(p1, p2, p4, p3, p5, p6, p7, lw=1, layout=l, size=(600, 800))
#savefig("freeze_activation_individual_self.tex")
savefig("jump_activation_individual_self.svg")

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

t₁ = 400.0
rawsim = simulate(valuemodel, Dict(:t₁ => t₁))
sim = link(valuemodel, rawsim)
ylabs = link_names(valuemodel)
param = param_value
#hillfun(t₁) = (x, t) -> 1 - hill(param[:kab], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
p1 = plot(times, sim[:,1], lab="", lc=:red, lw=2, yscale=:log10, ylab=ylabs[1])
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
p2 = plot(times, sim[:,2], lab="", lc=:red, lw=2, ylab=ylabs[2])
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
p3 = plot(times, param_rt[:wt], lab="wt", lc=:black)
plot!(p3, times, valuemodel_r(t₁).(times), lab="del", ylab="Activation rate", lc=:red)
p4 = plot(times, hillfun(800.0).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:black)
plot!(p4, times, hillfun(t₁).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:red)
plot(p1, p2, p4, p3, xlab="Time (days)", lw=2)
savefig("valuemodel_interaction.svg")

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

t₁ = 600.0
rawsim = simulate(valuemodel, Dict(:t₁ => t₁))
sim = link(valuemodel, rawsim)
ylabs = link_names(valuemodel)
param = param_value
#hillfun(t₁) = (x, t) -> 1 - hill(param[:kab], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
p1 = plot(times, sim[:,1], lab="", lc=:red, lw=2, yscale=:log10, ylab=ylabs[1])
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
p2 = plot(times, sim[:,2], lab="", lc=:red, lw=2, ylab=ylabs[2])
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
p3 = plot(times, param_rt[:wt], lab="wt", lc=:black)
plot!(p3, times, valuemodel_r(t₁).(times), lab="del", ylab="Activation rate", lc=:red)
p4 = plot(times, hillfun(800.0).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:black)
plot!(p4, times, hillfun(t₁).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:red)
plot(p1, p2, p4, p3, xlab="Time (days)", lw=2)

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

t₁ = 30.0
rawsim = simulate(valuemodel, Dict(:t₁ => t₁))
sim = link(valuemodel, rawsim)
ylabs = link_names(valuemodel)
param = param_value
#hillfun(t₁) = (x, t) -> 1 - hill(param[:kab], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
p1 = plot(times, sim[:,1], lab="", lc=:red, lw=2, yscale=:log10, ylab=ylabs[1])
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
p2 = plot(times, sim[:,2], lab="", lc=:red, lw=2, ylab=ylabs[2])
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
p3 = plot(times, param_rt[:wt], lab="wt", lc=:black)
plot!(p3, times, valuemodel_r(t₁).(times), lab="del", ylab="Activation rate", lc=:red)
p4 = plot(times, hillfun(800.0).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:black)
plot!(p4, times, hillfun(t₁).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:red)
plot(p1, p2, p4, p3, xlab="Time (days)", lw=2)

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

itimes = collect(0:10:699)#[100, 200]#, 60, 100, 300, 600]
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
rawsim = simulate(valuemodel, Dict(:t₁ => 700.0))
sim = link(valuemodel, rawsim)
ylabs = link_names(valuemodel)
p1 = plot(times, sim[:,1], lab="", lc=:gray, lw=2, yscale=:log10, ylab=ylabs[1])
p2 = plot(ylab=ylabs[2], xlab="Age (days)", ylim=(0.18, 0.65))
#p2 = plot(times, sim[:,2], lab="", lc=:gray, lw=2, ylab=ylabs[2])
p3 = plot(times, param_rt[:wt], lab="", ylab="Activation rate", lc=:gray, lw=2)
#p4 = plot(times, param_bt[:wt], lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:black, lw=2)
p4 = plot(times, hillfun(800).(sim[:,1], times), lab="", 
          ylab="Self-renewal probability", legend=:bottomright, lc=:black, lw=2)
for (i, t₁) in enumerate(itimes)
    rawsim = simulate(valuemodel, Dict(:t₁ => Float64(t₁)))
    sim = link(valuemodel, rawsim)
    has_changed = times .> t₁
    has_changed[max(1, findfirst(has_changed)-3):end] .= true
    plot!(p1, times[has_changed], sim[has_changed,1], lab="", lw=1, lc=colours[i]) #t₁=$t₁
    plot!(p2, times[has_changed], sim[has_changed,2], lab="", lw=1, lc=colours[i])
    plot!(p3, times[has_changed], valuemodel_r(t₁).(times[has_changed]), lab="", lw=1, lc=colours[i])
    plot!(p4, times[has_changed], hillfun(t₁).(sim[:,1], times)[has_changed], lab="", lw=1, lc=colours[i])
end
#plot!(p1, times, sim[:,1], lab="", lc=:red, lw=2)
@df @where(data, :genotype .== "wt", :name .== "total")  scatter!(p1, :age, :value, mc=:black, lab="")
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
plot(p1, p2, p4, p3, xlab="Time (days)", size=(600, 600))

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

itimes = collect(0:1:660)#[100, 200]#, 60, 100, 300, 600]
itimes = metrics_loss.itime
rawsim = simulate(valuemodel, Dict(:t₁ => 700.0))
regsim = link(valuemodel, rawsim)
nsc660 = regsim[times .== 660.0, 1][1]
nscs = map(itimes) do t₁
    rawsim = simulate(valuemodel, Dict(:t₁ => t₁))
    regsim = link(valuemodel, rawsim)
    nsc660 .- regsim[times .== 660.0, 1][1]
end
metrics_loss[!,"value_fullself"] = nscs
plot(itimes, nscs, lc=:black, lab="", xlab="Intervention Time (days)", ylab="Stem cells lost", title="Loss of stem cells from intervention")

itimes = collect(0:10:600)#[100, 200]#, 60, 100, 300, 600]
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
rawsim = simulate(valuemodel, Dict(:t₁ => 700.0))
sim = link(valuemodel, rawsim)
p = plot(xlab="Time (days)", ylab="Rate of TAP production from NSCs", yscale=:log10)
#hillfun(t₁) = (x, t) -> 1 - hill(param[:kab], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
for (i, t₁) in enumerate(itimes)
    rawsim = simulate(valuemodel, Dict(:t₁ => Float64(t₁)))
    #has_changed = times .> t₁
    plot!(p, times, rawsim[:,2] .* 2 .* 0.95 .* (1 .- hillfun(t₁).(rawsim[:,1], times)), lab="", lc=colours[i])
end
rawsim = simulate(valuemodel, Dict(:t₁ => Float64(700)))
plot!(p, times, rawsim[:,2] .* 2 .* 0.95 .* (1 .- param_bt[:wt].(times)), lab="WT", lc=:gray, lw=2)
#plot!(p, sim_ifnko_both.t, sim_ifnko_both.active .* sim_ifnko_both.counts .* sim_ifnko_both.b .* 2 .* 0.95, lab="KO", lc=:red, lw=2)
p
#savefig("plots/scenario1-lifelong-neurogenesis.svg")

itimes = metrics_loss.itime
int_neurogenesis_1 = zeros(Float64, length(itimes))
#int_neurogenesis_2 = zeros(Float64, length(itimes))
int_neurogenesis_wt = zeros(Float64, length(itimes))
step = times[2] - times[1]
#hillfun(t₁) = (x, t) -> 1 - hill(param[:kab], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
for (i, t₁) in enumerate(itimes)
    rawsim = simulate(valuemodel, Dict(:t₁ => Float64(t₁)))
    int_neurogenesis_1[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- hillfun(t₁).(rawsim[:,1], times)))
    #rawsim = simulate(valuemodel, Dict(:t₁ => Float64(t₁)))
    #int_neurogenesis_2[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- param_bt[:wt].(times)))
    rawsim = simulate(valuemodel, Dict(:t₁ => Float64(700)))
    int_neurogenesis_wt[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- hillfun(t₁).(rawsim[:,1], times)))
end
yticks = 50:5:115
metrics_neurogen[!,"value_fullself"] = int_neurogenesis_1 ./ int_neurogenesis_wt
plot(xlab="Intervention Time (days)", ylab="increase in TAPs produced over wt", title="Life-long (700 days) neurogenesis", yticks=(collect(yticks), ["$x\\%" for x in yticks]))
plot!(itimes, int_neurogenesis_1 ./ int_neurogenesis_wt .* 100, lab="'beta' model")
#plot!(itimes, int_neurogenesis_2 ./ int_neurogenesis_wt .* 100, lab="'value' model")
#plot!(itimes, int_neurogenesis_wt .* step, lab="wt", lc=:gray)

Self-renewal is partially genotype dependent

#=
param_value = Dict( # from cofit
    :nsc₀ => 1986.0,
    :r₀ => 0.8426,
    :r₁ => 0.4928,
    :βr => 0.001891,
    :kab => 141.987,
    :nb_wt => 0.0364225,
    :nb_ko => 0.0364225,
    #:nb_ko => 0.0765863,
    :pₛ => 0.95,
)
=#
param_value = Dict( # from cofit
    :nsc₀ => 3000.0,
    :r₀ => 0.839538,
    :r₁ => 0.491999,
    :βr => 0.00188605,
    :kab_wt => 153.81,
    :kab_ko => 153.81,
    :nb_wt => 0.039907,
    :nb_ko => 0.0841552,
    #:nb_ko => 0.0765863,
    :pₛ => 0.95,
)

Dict{Symbol, Float64} with 9 entries:
  :kab_ko => 153.81
  :kab_wt => 153.81
  :nb_ko  => 0.0841552
  :βr     => 0.00188605
  :r₀     => 0.839538
  :nsc₀   => 3000.0
  :r₁     => 0.491999
  :nb_wt  => 0.039907
  :pₛ     => 0.95

times = collect(0:1:700.0) # 0.1 was step
valuemodel = ChangedValueSelf((0.0, 700.0), times, zeros(Int64, length(times)), param_value)

ChangedValueSelf((0.0, 700.0), [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0  …  691.0, 692.0, 693.0, 694.0, 695.0, 696.0, 697.0, 698.0, 699.0, 700.0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0  …  0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [5 => 153.81, 4 => 153.81, 7 => 0.0841552, 3 => 0.00188605, 2 => 0.839538, 1 => 3000.0, 9 => 0.491999, 6 => 0.039907, 8 => 0.95])

free_parameters(valuemodel)

1-element Vector{Symbol}:
 :t₁

valuemodel_r = t₁ -> broken_exponential_decay(param_value[:r₀], param_value[:βr], param_value[:r₁], t₁)

#515 (generic function with 1 method)

sim

701×2 Matrix{Float64}:
 3000.0    0.918556
 2886.03   0.549385
 2805.52   0.480759
 2733.0    0.467792
 2663.92   0.464993
 2597.4    0.464051
 2533.19   0.463452
 2471.19   0.462916
 2411.29   0.462385
 2353.42   0.461863
 2297.48   0.461336
 2243.39   0.460822
 2191.09   0.460278
    ⋮      
  155.599  0.192851
  155.572  0.192567
  155.545  0.192274
  155.518  0.191983
  155.491  0.191729
  155.465  0.191459
  155.439  0.191129
  155.413  0.190802
  155.387  0.190644
  155.361  0.19043
  155.335  0.189941
  155.309  0.19011

#pgfplots()
pyplot()
itimes = collect(0:1:660)
rawsim = simulate(valuemodel, Dict(:t₁ => 800.0)) # No intervention model
nonsim = link(valuemodel, rawsim)
sim = nonsim
ylabs = link_names(valuemodel)
param = param_beta
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
# Plot base simulation (wt)
p1 = plot(yscale=:log10, ylab=ylabs[1], xlab="Age (days)")
p2 = plot(ylab="Relative NSCs", xlab="Age (days)")
p3 = plot(xlab="Age (days)", ylab="Activation rate")
p4 = plot(xlab="Age (days)",  ylab="Self-renewal probability", legend=:bottomright)
# ???
neurogen_wt = 2 .* 0.95 .* (1 .- hillfun(800.0).(sim[:,1], times)) .* rawsim[:,2]
tot = 2 .* 0.95 .* (1 .- hillfun(800.0).(sim[:,1], times)) .* rawsim[:,2]
tot_wt = sim[:,1]
tot_rel = DataFrame(:t => times, :wt => tot_wt ./ tot_wt)
neurogen_rel = DataFrame(:t => times, :wt => neurogen_wt ./ neurogen_wt)
tot_abs = DataFrame(:t => times, :wt => tot_wt)
neurogen_abs = DataFrame(:t => times, :wt => neurogen_wt)
p5 = plot(yscale=:log10, ylab="Rate of progenitor production", xlab="Age (days)")
p6 = plot(ylab="Relative rate", xlab="Age (days)")
# Plot interventions
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
# Compute some things
step = times[2] - times[1]
non_neurogen = step .* neurogen_wt
# Preallocate vectors
loss = zeros(Float64, length(itimes))
total_neurogen = zeros(Float64, length(itimes))
for (i, itime) in enumerate(itimes)
    rawsim = simulate(valuemodel, Dict(:t₁ => Float64(itime)))
    changed = times .> itime
    sim = link(valuemodel, rawsim)
    loss[i] = nonsim[findfirst(times .== 660.0), 1] .- sim[findfirst(times .== 660.0),1]
    neurogen = 2 .* 0.95 .* (1 .- hillfun(itime).(sim[:,1], times)) .* rawsim[:,2]
    total_neurogen[i] = sum(neurogen .* step) ./ sum(non_neurogen)
    plot!(p1, times[changed], sim[changed,1], lab="", lc=colours[i], la=0.9)
    plot!(p2, times, sim[:,1] ./ tot_wt, lab="", lc=colours[i], la=0.9)
    #plot!(p3, times[changed], rawsim[changed,3], lab="", lc=colours[i], la=0.9)
    plot!(p3, times[changed], valuemodel_r(itime).(times[changed]), lab="", lc=colours[i], la=0.9)
    plot!(p4, times[changed], hillfun(itime).(sim[changed,1], times[changed]), lab="", lc=colours[i], la=0.9)
    plot!(p5, times[changed], neurogen[changed], lab="", lc=colours[i], la=0.9)
    plot!(p6, times, neurogen ./ neurogen_wt, lab="", lc=colours[i], la=0.9)
    tot_rel[!,"d$(Int64(round(itime)))"] = sim[:,1] ./ tot_wt
    neurogen_rel[!,"d$(Int64(round(itime)))"] = neurogen ./ neurogen_wt
    tot_abs[!,"d$(Int64(round(itime)))"] = sim[:,1]
    neurogen_abs[!,"d$(Int64(round(itime)))"] = neurogen
end
# Plot previous simulation
plot!(p1, times, nonsim[:,1], lab="", lc=:black, lw=2)
plot!(p2, times, nonsim[:,1] ./ nonsim[:,1], lab="", lc=:black, lw=2)
plot!(p3, times, valuemodel_r(900.0), lab="", lc=:black, lw=2)
plot!(p4, times, hillfun(800.0).(nonsim[:,1], times), lab="", lc=:black, lw=2)
plot!(p5, times, neurogen_wt, lc=:black, lw=2, lab="")
plot!(p6, times, neurogen_wt ./ neurogen_wt, lc=:black, lw=2, lab="")
# Plot data
#@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
#@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
# Plot metrics
p7 = plot(itimes, loss, lab="", ylab="Stem cells lost", xlab="Intervention time (days)", lc=:black)
steps = round(minimum(total_neurogen); sigdigits=1):0.05:round(maximum(total_neurogen); sigdigits=1)
p8 = plot(itimes, total_neurogen, lc=:black, lab="", ylab="Total progenitor gain", xlab = "Intervention time (days)",
          ytick=(collect(steps), ["$(Int64(100*x))\\%" for x in steps]))
#p7 = plot(itimes, total_neurogen)
l = @layout [
    p1 p2;
    p3 p4;
    p5 p6;
    p7 p8;
]
plot(p1, p2, p4, p3, p5, p6, p7, p8, lw=1, layout=l, size=(600, 800))
#savefig("jump_activation_semishared_self.tex")
savefig("jump_activation_semishared_self.svg")
#tot_rel |> CSV.write("sims/jump_activation_semishared_self_relative_nsc.csv")
#neurogen_rel |> CSV.write("sims/jump_activation_semishared_self_relative_neurogen.csv")
#tot_abs |> CSV.write("sims/jump_activation_semishared_self_absolute_nsc.csv")
#neurogen_abs |> CSV.write("sims/jump_activation_semishared_self_absolute_neurogen.csv")

┌ Warning: No strict ticks found
└ @ PlotUtils /home/jooa/.julia/packages/PlotUtils/VgXdq/src/ticks.jl:294

#pgfplots()
pyplot()
rawsim = simulate(valuemodel, Dict(:t₁ => 800.0)) # No intervention model
itimes = collect(0:10:660)
nonsim = link(valuemodel, rawsim)
ylabs = link_names(valuemodel)
param = param_beta
hillfun(t₁) = (x, t) -> 1 - hill(param[:kab_wt], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
# Plot base simulation (wt)
p1 = plot(yscale=:log10, ylab=ylabs[1], xlab="Age (days)")
#p2 = plot(ylab=ylabs[2], xlab="Age (days)")
p2 = plot(ylab=ylabs[2], xlab="Age (days)", ylim=(0.18, 0.65))
p3 = plot(xlab="Age (days)", ylab="Activation rate")
p4 = plot(xlab="Age (days)",  ylab="Self-renewal probability", legend=:bottomright)
# ???
neurogen = 2 .* 0.95 .* (1 .- hillfun(800.0).(sim[:,1], times)) .* rawsim[:,2]
p5 = plot(yscale=:log10, ylab="Rate of progenitor production", xlab="Age (days)")
# Plot interventions
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
# Compute some things
step = times[2] - times[1]
non_neurogen = step .* neurogen
# Preallocate vectors
loss = zeros(Float64, length(itimes))
total_neurogen = zeros(Float64, length(itimes))
for (i, itime) in enumerate(itimes)
    rawsim = simulate(valuemodel, Dict(:t₁ => Float64(itime)))
    changed = times .> itime
    changed[max(1, findfirst(changed)-5):end] .= true
    sim = link(valuemodel, rawsim)
    loss[i] = nonsim[findfirst(times .== 660.0), 1] .- sim[findfirst(times .== 660.0),1]
    neurogen_tmp = 2 .* 0.95 .* (1 .- hillfun(itime).(sim[:,1], times)) .* rawsim[:,2]
    total_neurogen[i] = sum(neurogen_tmp .* step) ./ sum(non_neurogen)
    plot!(p1, times[changed], sim[changed,1], lab="", lc=colours[i], la=0.9)
    plot!(p2, times[changed], sim[changed,2], lab="", lc=colours[i], la=0.9)
#    plot!(p3, times[changed], rawsim[changed,3], lab="", lc=colours[i], la=0.9)
    plot!(p3, times[changed], valuemodel_r(itime).(times[changed]), lab="", lc=colours[i], la=0.9)
    plot!(p4, times[changed], hillfun(itime).(sim[changed,1], times[changed]), lab="", lc=colours[i], la=0.9)
    plot!(p5, times[changed], neurogen_tmp[changed], lab="", lc=colours[i], la=0.9)
end
# Plot previous sims
plot!(p1, times, nonsim[:,1], lab="", lc=:black, lw=2)
plot!(p2, times, nonsim[:,2], lab="", lc=:black, lw=2)
#p3 = plot(times, param_rt[:wt], lab="wt", lc=:black)
#plot!(p3, times, valuemodel_r(t₁).(times), lab="del", ylab="Activation rate", lc=:red)
#plot!(p3, times, param_rt[:wt], lab="", lc=:black, lw=2)
plot!(p3, times, valuemodel_r(900.0), lab="", lc=:black, lw=2)
plot!(p4, times, hillfun(800.0).(nonsim[:,1], times), lab="", lc=:black, lw=2)
plot!(p5, times, neurogen, lc=:black, lw=2,lab="")
# Plot data
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
# Plot metrics
p6 = plot(itimes, loss, lab="", ylab="Stem cells lost", xlab="Intervention time (days)", lc=:black)
neurogen_range = maximum(total_neurogen) - minimum(total_neurogen)
steps = 0.80:0.05:1.1
p7 = plot(itimes, total_neurogen, lc=:black, lab="", ylab="Total progenitor gain", xlab = "Intervention time (days)",
          ytick=(collect(steps), ["$(Int64(round(100*x)))\\%" for x in steps]), ylim=(minimum(steps), maximum(steps)))
#p7 = plot(itimes, total_neurogen)
l = @layout [
    p1 p2;
    p3 p4;
    p5{.5h} [p6;
             p7]
]
plot(p1, p2, p4, p3, p5, p6, p7, lw=1, layout=l, size=(600, 800))
#savefig("freeze_activation_individual_self.tex")
savefig("jump_activation_semishared_self.svg")

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

total_mat = 0:10:660
#pgfplots()
pyplot()
rawsim = simulate(valuemodel, Dict(:t₁ => 800.0)) # No intervention model
itimes = collect(0:10:660)
all_nsc = zeros(length(times), length(itimes))
all_neurogen = zeros(length(times), length(itimes))
nonsim = link(valuemodel, rawsim)
ylabs = link_names(valuemodel)
param = param_beta
hillfun(t₁) = (x, t) -> 1 - hill(param[:kab_wt], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
# Plot base simulation (wt)
p1 = plot(yscale=:log10, ylab=ylabs[1], xlab="Age (days)")
#p2 = plot(ylab=ylabs[2], xlab="Age (days)")
p2 = plot(ylab=ylabs[2], xlab="Age (days)", ylim=(0.18, 0.65))
p3 = plot(xlab="Age (days)", ylab="Activation rate")
p4 = plot(xlab="Age (days)",  ylab="Self-renewal probability", legend=:bottomright)
# ???
neurogen = 2 .* 0.95 .* (1 .- hillfun(800.0).(sim[:,1], times)) .* rawsim[:,2]
p5 = plot(yscale=:log10, ylab="Rate of progenitor production", xlab="Age (days)")
# Plot interventions
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
# Compute some things
step = times[2] - times[1]
non_neurogen = step .* neurogen
# Preallocate vectors
loss = zeros(Float64, length(itimes))
total_neurogen = zeros(Float64, length(itimes))
for (i, itime) in enumerate(itimes)
    rawsim = simulate(valuemodel, Dict(:t₁ => Float64(itime)))
    changed = times .> itime
    changed[max(1, findfirst(changed)-5):end] .= true
    sim = link(valuemodel, rawsim)
    loss[i] = nonsim[findfirst(times .== 660.0), 1] .- sim[findfirst(times .== 660.0),1]
    neurogen_tmp = 2 .* 0.95 .* (1 .- hillfun(itime).(sim[:,1], times)) .* rawsim[:,2]
    total_neurogen[i] = sum(neurogen_tmp .* step) ./ sum(non_neurogen)
    all_nsc[:,i] = sim[:,1]
    plot!(p1, times[changed], sim[changed,1], lab="", lc=colours[i], la=0.9)
    plot!(p2, times[changed], sim[changed,2], lab="", lc=colours[i], la=0.9)
#    plot!(p3, times[changed], rawsim[changed,3], lab="", lc=colours[i], la=0.9)
    plot!(p3, times[changed], valuemodel_r(itime).(times[changed]), lab="", lc=colours[i], la=0.9)
    plot!(p4, times[changed], hillfun(itime).(sim[changed,1], times[changed]), lab="", lc=colours[i], la=0.9)
    plot!(p5, times[changed], neurogen_tmp[changed], lab="", lc=colours[i], la=0.9)
    all_neurogen[:,i] = neurogen_tmp
end
# Plot previous sims
plot!(p1, times, nonsim[:,1], lab="", lc=:black, lw=2)
plot!(p2, times, nonsim[:,2], lab="", lc=:black, lw=2)
#p3 = plot(times, param_rt[:wt], lab="wt", lc=:black)
#plot!(p3, times, valuemodel_r(t₁).(times), lab="del", ylab="Activation rate", lc=:red)
#plot!(p3, times, param_rt[:wt], lab="", lc=:black, lw=2)
plot!(p3, times, valuemodel_r(900.0), lab="", lc=:black, lw=2)
plot!(p4, times, hillfun(800.0).(nonsim[:,1], times), lab="", lc=:black, lw=2)
plot!(p5, times, neurogen, lc=:black, lw=2,lab="")
# Plot data
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
# Plot metrics
p6 = plot(itimes, loss, lab="", ylab="Stem cells lost", xlab="Intervention time (days)", lc=:black)
neurogen_range = maximum(total_neurogen) - minimum(total_neurogen)
steps = 0.80:0.05:1.1
p7 = plot(itimes, total_neurogen, lc=:black, lab="", ylab="Total progenitor gain", xlab = "Intervention time (days)",
          ytick=(collect(steps), ["$(Int64(round(100*x)))\\%" for x in steps]), ylim=(minimum(steps), maximum(steps)))
#p7 = plot(itimes, total_neurogen)
l = @layout [
    p1 p2;
    p3 p4;
    p5{.5h} [p6;
             p7]
]
plot(p1, p2, p4, p3, p5, p6, p7, lw=1, layout=l, size=(600, 800))
#savefig("freeze_activation_individual_self.tex")
savefig("jump_activation_semishared_self.svg")

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

DataFrame("time" => itimes, "loss" => loss, "neurogen" => total_neurogen) |> CSV.write("data_export/all_metrics.csv")

"data_export/all_metrics.csv"

neurogen

701-element Vector{Float64}:
 2905.8763736291944
 1661.6974622621224
 1416.2175028933223
 1343.0954577345854
 1300.6694641572635
 1264.6248548525637
 1230.7124056893203
 1198.1653243774574
 1166.8018276690952
 1136.522877208645
 1107.3733095659613
 1079.129443928931
 1052.0172627494176
    ⋮
   29.699723058277534
   29.652592643169385
   29.646707073288386
   29.560277645023245
   29.487034297870398
   29.48324653841649
   29.407130243155517
   29.350600296197122
   29.32605394129574
   29.254411321276336
   29.20846653849023
   29.16416832982201

all_neurogen_df = DataFrame(all_neurogen, ["d$i" for i in itimes])
all_neurogen_df[!,"wt"] = neurogen
all_neurogen_df[!,"time"] = times

701-element Vector{Float64}:
   0.0
   1.0
   2.0
   3.0
   4.0
   5.0
   6.0
   7.0
   8.0
   9.0
  10.0
  11.0
  12.0
   ⋮
 689.0
 690.0
 691.0
 692.0
 693.0
 694.0
 695.0
 696.0
 697.0
 698.0
 699.0
 700.0

all_nsc_df = DataFrame(all_nsc, ["d$i" for i in itimes])
all_nsc_df[!,"wt"] = nonsim[:,1] 
all_nsc_df[!,"time"] = times

701-element Vector{Float64}:
   0.0
   1.0
   2.0
   3.0
   4.0
   5.0
   6.0
   7.0
   8.0
   9.0
  10.0
  11.0
  12.0
   ⋮
 689.0
 690.0
 691.0
 692.0
 693.0
 694.0
 695.0
 696.0
 697.0
 698.0
 699.0
 700.0

all_neurogen_df |> CSV.write("data_export/all_neurogen.csv")
all_nsc_df |> CSV.write("data_export/all_nsc.csv")

"data_export/all_nsc.csv"

heatmap(log10.(all_nsc))

plot(neurogen, yscale=:log10)

itimes = collect(0:1:660)#[100, 200]#, 60, 100, 300, 600]
itimes = metrics_loss.itime
rawsim = simulate(valuemodel, Dict(:t₁ => 700.0))
regsim = link(valuemodel, rawsim)
nsc660 = regsim[times .== 660.0, 1][1]
nscs = map(itimes) do t₁
    rawsim = simulate(valuemodel, Dict(:t₁ => t₁))
    regsim = link(valuemodel, rawsim)
    nsc660 .- regsim[times .== 660.0, 1][1]
end
metrics_loss[!,"value_partself"] = nscs
plot(itimes, nscs, lc=:black, lab="", xlab="Intervention Time (days)", ylab="Stem cells lost", title="Loss of stem cells from intervention")

itimes = metrics_loss.itime
int_neurogenesis_1 = zeros(Float64, length(itimes))
#int_neurogenesis_2 = zeros(Float64, length(itimes))
int_neurogenesis_wt = zeros(Float64, length(itimes))
step = times[2] - times[1]
hillfun(t₁) = (x, t) -> 1 - hill(param[:kab_wt], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
for (i, t₁) in enumerate(itimes)
    rawsim = simulate(valuemodel, Dict(:t₁ => Float64(t₁)))
    int_neurogenesis_1[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- hillfun(t₁).(rawsim[:,1], times)))
    #rawsim = simulate(valuemodel, Dict(:t₁ => Float64(t₁)))
    #int_neurogenesis_2[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- param_bt[:wt].(times)))
    rawsim = simulate(valuemodel, Dict(:t₁ => Float64(700)))
    int_neurogenesis_wt[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- hillfun(t₁).(rawsim[:,1], times)))
end
metrics_neurogen[!,"value_partself"] = int_neurogenesis_1 ./ int_neurogenesis_wt
plot(xlab="Intervention Time (days)", ylab="increase in TAPs produced over wt", title="Life-long (700 days) neurogenesis", yticks=([100, 105, 110, 115], ["100\\%", "105\\%", "110\\%", "115\\%"]))
plot!(itimes, int_neurogenesis_1 ./ int_neurogenesis_wt .* 100, lab="'beta' model")
#plot!(itimes, int_neurogenesis_2 ./ int_neurogenesis_wt .* 100, lab="'value' model")
#plot!(itimes, int_neurogenesis_wt .* step, lab="wt", lc=:gray)

Self-renewal is genotype-invariant

#=
param_value = Dict( # from cofit
    :nsc₀ => 1986.0,
    :r₀ => 0.8426,
    :r₁ => 0.4928,
    :βr => 0.001891,
    :kab => 141.987,
    :nb_wt => 0.0364225,
    :nb_ko => 0.0364225,
    #:nb_ko => 0.0765863,
    :pₛ => 0.95,
)
=#
param_value = Dict( # from cofit
    :nsc₀ => 5000.0,
    :r₀ => 0.672387,
    :r₁ => 0.495324,
    :βr => 0.0015538,
    :kab_wt => 153.724,
    :kab_ko => 153.724,
    :nb_wt => 0.0584954,
    :nb_ko => 0.0584954,
    #:nb_ko => 0.0765863,
    :pₛ => 0.95,
)

Dict{Symbol, Float64} with 9 entries:
  :kab_ko => 153.724
  :kab_wt => 153.724
  :nb_ko  => 0.0584954
  :βr     => 0.0015538
  :r₀     => 0.672387
  :nsc₀   => 5000.0
  :r₁     => 0.495324
  :nb_wt  => 0.0584954
  :pₛ     => 0.95

times = collect(0:1:700.0)
valuemodel = ChangedValueSelf((0.0, 700.0), times, zeros(Int64, length(times)), param_value)

ChangedValueSelf((0.0, 700.0), [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0  …  691.0, 692.0, 693.0, 694.0, 695.0, 696.0, 697.0, 698.0, 699.0, 700.0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0  …  0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [5 => 153.724, 4 => 153.724, 7 => 0.0584954, 3 => 0.0015538, 2 => 0.672387, 1 => 5000.0, 9 => 0.495324, 6 => 0.0584954, 8 => 0.95])

free_parameters(valuemodel)

1-element Vector{Symbol}:
 :t₁

valuemodel_r = t₁ -> broken_exponential_decay(param_value[:r₀], param_value[:βr], param_value[:r₁], t₁)

#561 (generic function with 1 method)

#pgfplots()
pyplot()
itimes = collect(0:1:660)
rawsim = simulate(valuemodel, Dict(:t₁ => 800.0)) # No intervention model
nonsim = link(valuemodel, rawsim)
sim = nonsim
ylabs = link_names(valuemodel)
param = param_beta
hillfun(t₁) = (x, t) -> 1 - hill(t < t₁ ? param[:kab_wt] : param[:kab_ko], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
# Plot base simulation (wt)
p1 = plot(yscale=:log10, ylab=ylabs[1], xlab="Age (days)")
p2 = plot(ylab="Relative NSCs", xlab="Age (days)")
p3 = plot(xlab="Age (days)", ylab="Activation rate")
p4 = plot(xlab="Age (days)",  ylab="Self-renewal probability", legend=:bottomright)
# ???
neurogen_wt = 2 .* 0.95 .* (1 .- hillfun(800.0).(sim[:,1], times)) .* rawsim[:,2]
tot = 2 .* 0.95 .* (1 .- hillfun(800.0).(sim[:,1], times)) .* rawsim[:,2]
tot_wt = sim[:,1]
tot_rel = DataFrame(:t => times, :wt => tot_wt ./ tot_wt)
neurogen_rel = DataFrame(:t => times, :wt => neurogen_wt ./ neurogen_wt)
tot_abs = DataFrame(:t => times, :wt => tot_wt)
neurogen_abs = DataFrame(:t => times, :wt => neurogen_wt)
p5 = plot(yscale=:log10, ylab="Rate of progenitor production", xlab="Age (days)")
p6 = plot(ylab="Relative rate", xlab="Age (days)")
# Plot interventions
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
# Compute some things
step = times[2] - times[1]
non_neurogen = step .* neurogen_wt
# Preallocate vectors
loss = zeros(Float64, length(itimes))
total_neurogen = zeros(Float64, length(itimes))
for (i, itime) in enumerate(itimes)
    rawsim = simulate(valuemodel, Dict(:t₁ => Float64(itime)))
    changed = times .> itime
    sim = link(valuemodel, rawsim)
    loss[i] = nonsim[findfirst(times .== 660.0), 1] .- sim[findfirst(times .== 660.0),1]
    neurogen = 2 .* 0.95 .* (1 .- hillfun(itime).(sim[:,1], times)) .* rawsim[:,2]
    total_neurogen[i] = sum(neurogen .* step) ./ sum(non_neurogen)
    plot!(p1, times[changed], sim[changed,1], lab="", lc=colours[i], la=0.9)
    plot!(p2, times, sim[:,1] ./ tot_wt, lab="", lc=colours[i], la=0.9)
    #plot!(p3, times[changed], rawsim[changed,3], lab="", lc=colours[i], la=0.9)
    plot!(p3, times[changed], valuemodel_r(itime).(times[changed]), lab="", lc=colours[i], la=0.9)
    plot!(p4, times[changed], hillfun(itime).(sim[changed,1], times[changed]), lab="", lc=colours[i], la=0.9)
    plot!(p5, times[changed], neurogen[changed], lab="", lc=colours[i], la=0.9)
    plot!(p6, times, neurogen ./ neurogen_wt, lab="", lc=colours[i], la=0.9)
    tot_rel[!,"d$(Int64(round(itime)))"] = sim[:,1] ./ tot_wt
    neurogen_rel[!,"d$(Int64(round(itime)))"] = neurogen ./ neurogen_wt
    tot_abs[!,"d$(Int64(round(itime)))"] = sim[:,1]
    neurogen_abs[!,"d$(Int64(round(itime)))"] = neurogen
end
# Plot previous simulation
plot!(p1, times, nonsim[:,1], lab="", lc=:black, lw=2)
plot!(p2, times, nonsim[:,1] ./ nonsim[:,1], lab="", lc=:black, lw=2)
plot!(p3, times, valuemodel_r(900.0), lab="", lc=:black, lw=2)
plot!(p4, times, hillfun(800.0).(nonsim[:,1], times), lab="", lc=:black, lw=2)
plot!(p5, times, neurogen_wt, lc=:black, lw=2, lab="")
plot!(p6, times, neurogen_wt ./ neurogen_wt, lc=:black, lw=2, lab="")
# Plot data
#@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
#@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
# Plot metrics
p7 = plot(itimes, loss, lab="", ylab="Stem cells lost", xlab="Intervention time (days)", lc=:black)
steps = round(minimum(total_neurogen); sigdigits=1):0.05:round(maximum(total_neurogen); sigdigits=1)
p8 = plot(itimes, total_neurogen, lc=:black, lab="", ylab="Total progenitor gain", xlab = "Intervention time (days)",
          ytick=(collect(steps), ["$(Int64(100*x))\\%" for x in steps]))
#p7 = plot(itimes, total_neurogen)
l = @layout [
    p1 p2;
    p3 p4;
    p5 p6;
    p7 p8;
]
plot(p1, p2, p4, p3, p5, p6, p7, p8, lw=1, layout=l, size=(600, 800))
#savefig("jump_activation_shared_self.tex")
savefig("jump_activation_shared_self.svg")
tot_rel |> CSV.write("sims/jump_activation_shared_self_relative_nsc.csv")
neurogen_rel |> CSV.write("sims/jump_activation_shared_self_relative_neurogen.csv")
tot_abs |> CSV.write("sims/jump_activation_shared_self_absolute_nsc.csv")
neurogen_abs |> CSV.write("sims/jump_activation_shared_self_absolute_neurogen.csv")

┌ Warning: No strict ticks found
└ @ PlotUtils /home/jooa/.julia/packages/PlotUtils/VgXdq/src/ticks.jl:294

"sims/jump_activation_shared_self_absolute_neurogen.csv"

#pgfplots()
pyplot()
rawsim = simulate(valuemodel, Dict(:t₁ => 800.0)) # No intervention model
itimes = collect(0:10:660)
nonsim = link(valuemodel, rawsim)
ylabs = link_names(valuemodel)
param = param_beta
hillfun(t₁) = (x, t) -> 1 - hill(param[:kab_wt], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
# Plot base simulation (wt)
p1 = plot(yscale=:log10, ylab=ylabs[1], xlab="Age (days)")
#p2 = plot(ylab=ylabs[2], xlab="Age (days)")
p2 = plot(ylab=ylabs[2], xlab="Age (days)", ylim=(0.18, 0.65))
p3 = plot(xlab="Age (days)", ylab="Activation rate")
p4 = plot(xlab="Age (days)",  ylab="Self-renewal probability", legend=:bottomright)
# ???
neurogen = 2 .* 0.95 .* (1 .- hillfun(800.0).(sim[:,1], times)) .* rawsim[:,2]
p5 = plot(yscale=:log10, ylab="Rate of progenitor production", xlab="Age (days)")
# Plot interventions
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
# Compute some things
step = times[2] - times[1]
non_neurogen = step .* neurogen
# Preallocate vectors
loss = zeros(Float64, length(itimes))
total_neurogen = zeros(Float64, length(itimes))
for (i, itime) in enumerate(itimes)
    rawsim = simulate(valuemodel, Dict(:t₁ => Float64(itime)))
    changed = times .> itime
    changed[max(1, findfirst(changed)-5):end] .= true
    sim = link(valuemodel, rawsim)
    loss[i] = nonsim[findfirst(times .== 660.0), 1] .- sim[findfirst(times .== 660.0),1]
    neurogen = 2 .* 0.95 .* (1 .- hillfun(itime).(sim[:,1], times)) .* rawsim[:,2]
    total_neurogen[i] = sum(neurogen .* step) ./ sum(non_neurogen)
    plot!(p1, times[changed], sim[changed,1], lab="", lc=colours[i], la=0.9)
    plot!(p2, times[changed], sim[changed,2], lab="", lc=colours[i], la=0.9)
#    plot!(p3, times[changed], rawsim[changed,3], lab="", lc=colours[i], la=0.9)
    plot!(p3, times[changed], valuemodel_r(itime).(times[changed]), lab="", lc=colours[i], la=0.9)
    plot!(p4, times[changed], hillfun(itime).(sim[changed,1], times[changed]), lab="", lc=colours[i], la=0.9)
    plot!(p5, times[changed], neurogen[changed], lab="", lc=colours[i], la=0.9)
end
# Plot previous sims
plot!(p1, times, nonsim[:,1], lab="", lc=:black, lw=2)
plot!(p2, times, nonsim[:,2], lab="", lc=:black, lw=2)
#p3 = plot(times, param_rt[:wt], lab="wt", lc=:black)
#plot!(p3, times, valuemodel_r(t₁).(times), lab="del", ylab="Activation rate", lc=:red)
#plot!(p3, times, param_rt[:wt], lab="", lc=:black, lw=2)
plot!(p3, times, valuemodel_r(900.0), lab="", lc=:black, lw=2)
plot!(p4, times, hillfun(800.0).(nonsim[:,1], times), lab="", lc=:black, lw=2)
plot!(p5, times, neurogen, lc=:black, lw=2,lab="")
# Plot data
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
# Plot metrics
p6 = plot(itimes, loss, lab="", ylab="Stem cells lost", xlab="Intervention time (days)", lc=:black)
neurogen_range = maximum(total_neurogen) - minimum(total_neurogen)
steps = 1.00:0.05:1.1
p7 = plot(itimes, total_neurogen, lc=:black, lab="", ylab="Total progenitor gain", xlab = "Intervention time (days)",
          ytick=(collect(steps), ["$(Int64(round(100*x)))\\%" for x in steps]))
#p7 = plot(itimes, total_neurogen)
l = @layout [
    p1 p2;
    p3 p4;
    p5{.5h} [p6;
             p7]
]
plot(p1, p2, p4, p3, p5, p6, p7, lw=1, layout=l, size=(600, 800))
#savefig("freeze_activation_individual_self.tex")
savefig("jump_activation_shared_self.svg")

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

t₁ = 400.0
rawsim = simulate(valuemodel, Dict(:t₁ => t₁))
sim = link(valuemodel, rawsim)
ylabs = link_names(valuemodel)
param = param_value
hillfun(t₁) = (x, t) -> 1 - hill(param[:kab_wt], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
p1 = plot(times, sim[:,1], lab="", lc=:red, lw=2, yscale=:log10, ylab=ylabs[1])
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
p2 = plot(times, sim[:,2], lab="", lc=:red, lw=2, ylab=ylabs[2])
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
p3 = plot(times, param_rt[:wt], lab="wt", lc=:black)
plot!(p3, times, valuemodel_r(t₁).(times), lab="del", ylab="Activation rate", lc=:red)
p4 = plot(times, hillfun(800.0).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:black)
plot!(p4, times, hillfun(t₁).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:red)
plot(p1, p2, p4, p3, xlab="Time (days)", lw=2)
savefig("valuemodel_static.svg")

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

t₁ = 600.0
rawsim = simulate(valuemodel, Dict(:t₁ => t₁))
sim = link(valuemodel, rawsim)
ylabs = link_names(valuemodel)
param = param_value
hillfun(t₁) = (x, t) -> 1 - hill(param[:kab_wt], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
p1 = plot(times, sim[:,1], lab="", lc=:red, lw=2, yscale=:log10, ylab=ylabs[1])
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
p2 = plot(times, sim[:,2], lab="", lc=:red, lw=2, ylab=ylabs[2])
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
p3 = plot(times, param_rt[:wt], lab="wt", lc=:black)
plot!(p3, times, valuemodel_r(t₁).(times), lab="del", ylab="Activation rate", lc=:red)
p4 = plot(times, hillfun(800.0).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:black)
plot!(p4, times, hillfun(t₁).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:red)
plot(p1, p2, p4, p3, xlab="Time (days)", lw=2)

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

t₁ = 30.0
rawsim = simulate(valuemodel, Dict(:t₁ => t₁))
sim = link(valuemodel, rawsim)
ylabs = link_names(valuemodel)
param = param_value
hillfun(t₁) = (x, t) -> 1 - hill(param[:kab_wt], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
p1 = plot(times, sim[:,1], lab="", lc=:red, lw=2, yscale=:log10, ylab=ylabs[1])
@df @where(data, :genotype .== "wt", :name .== "total") scatter!(p1, :age, :value, mc=:black, lab="")
p2 = plot(times, sim[:,2], lab="", lc=:red, lw=2, ylab=ylabs[2])
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
p3 = plot(times, param_rt[:wt], lab="wt", lc=:black)
plot!(p3, times, valuemodel_r(t₁).(times), lab="del", ylab="Activation rate", lc=:red)
p4 = plot(times, hillfun(800.0).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:black)
plot!(p4, times, hillfun(t₁).(sim[:,1], times), lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:red)
plot(p1, p2, p4, p3, xlab="Time (days)", lw=2)

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

itimes = collect(0:10:699)#[100, 200]#, 60, 100, 300, 600]
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
rawsim = simulate(valuemodel, Dict(:t₁ => 700.0))
sim = link(valuemodel, rawsim)
ylabs = link_names(valuemodel)
p1 = plot(times, sim[:,1], lab="", lc=:gray, lw=2, yscale=:log10, ylab=ylabs[1])
p2 = plot(times, sim[:,2], lab="", lc=:gray, lw=2, ylab=ylabs[2])
p3 = plot(times, param_rt[:wt], lab="", ylab="Activation rate", lc=:gray, lw=2)
#p4 = plot(times, param_bt[:wt], lab="", ylab="Self-renewal probability", legend=:bottomright, lc=:black, lw=2)
p4 = plot(times, hillfun(800).(sim[:,1], times), lab="", 
          ylab="Self-renewal probability", legend=:bottomright, lc=:black, lw=2)
for (i, t₁) in enumerate(itimes)
    rawsim = simulate(valuemodel, Dict(:t₁ => Float64(t₁)))
    sim = link(valuemodel, rawsim)
    has_changed = times .> t₁
    has_changed[max(1, findfirst(has_changed)-3):end] .= true
    plot!(p1, times[has_changed], sim[has_changed,1], lab="", lw=1, lc=colours[i]) #t₁=$t₁
    plot!(p2, times[has_changed], sim[has_changed,2], lab="", lw=1, lc=colours[i])
    plot!(p3, times[has_changed], valuemodel_r(t₁).(times[has_changed]), lab="", lw=1, lc=colours[i])
    plot!(p4, times[has_changed], hillfun(t₁).(sim[:,1], times)[has_changed], lab="", lw=1, lc=colours[i])
end
#plot!(p1, times, sim[:,1], lab="", lc=:red, lw=2)
@df @where(data, :genotype .== "wt", :name .== "total")  scatter!(p1, :age, :value, mc=:black, lab="")
@df @where(data, :genotype .== "wt", :name .== "active") scatter!(p2, :age, :value, mc=:black, lab="")
plot(p1, p2, p4, p3, xlab="Time (days)", size=(600, 600))

┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822
┌ Warning: `@where is deprecated, use `@subset` instead.
└ @ DataFramesMeta /home/jooa/.julia/packages/DataFramesMeta/yzaoq/src/macros.jl:822

itimes = collect(0:1:660)#[100, 200]#, 60, 100, 300, 600]
itimes = metrics_loss.itime
rawsim = simulate(valuemodel, Dict(:t₁ => 700.0))
regsim = link(valuemodel, rawsim)
nsc660 = regsim[times .== 660.0, 1][1]
nscs = map(itimes) do t₁
    rawsim = simulate(valuemodel, Dict(:t₁ => t₁))
    regsim = link(valuemodel, rawsim)
    nsc660 .- regsim[times .== 660.0, 1][1]
end
metrics_loss[!,"value_noself"] = nscs
plot(itimes, nscs, lc=:black, lab="", xlab="Intervention Time (days)", ylab="Stem cells lost", title="Loss of stem cells from intervention")

itimes = collect(0:10:600)#[100, 200]#, 60, 100, 300, 600]
colours = cgrad(:roma, length(itimes), categorical=true, rev=false)
rawsim = simulate(valuemodel, Dict(:t₁ => 700.0))
sim = link(valuemodel, rawsim)
p = plot(xlab="Time (days)", ylab="Rate of TAP production from NSCs", yscale=:log10)
hillfun(t₁) = (x, t) -> 1 - hill(param[:kab_wt], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
for (i, t₁) in enumerate(itimes)
    rawsim = simulate(valuemodel, Dict(:t₁ => Float64(t₁)))
    #has_changed = times .> t₁
    plot!(p, times, rawsim[:,2] .* 2 .* 0.95 .* (1 .- hillfun(t₁).(rawsim[:,1], times)), lab="", lc=colours[i])
end
rawsim = simulate(valuemodel, Dict(:t₁ => Float64(700)))
plot!(p, times, rawsim[:,2] .* 2 .* 0.95 .* (1 .- param_bt[:wt].(times)), lab="WT", lc=:gray, lw=2)
#plot!(p, sim_ifnko_both.t, sim_ifnko_both.active .* sim_ifnko_both.counts .* sim_ifnko_both.b .* 2 .* 0.95, lab="KO", lc=:red, lw=2)
p
#savefig("plots/scenario1-lifelong-neurogenesis.svg")

itimes = metrics_loss.itime
int_neurogenesis_1 = zeros(Float64, length(itimes))
#int_neurogenesis_2 = zeros(Float64, length(itimes))
int_neurogenesis_wt = zeros(Float64, length(itimes))
step = times[2] - times[1]
hillfun(t₁) = (x, t) -> 1 - hill(param[:kab_wt], t < t₁ ? param[:nb_wt] : param[:nb_ko])(x)
for (i, t₁) in enumerate(itimes)
    rawsim = simulate(valuemodel, Dict(:t₁ => Float64(t₁)))
    int_neurogenesis_1[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- hillfun(t₁).(rawsim[:,1], times)))
    #rawsim = simulate(valuemodel, Dict(:t₁ => Float64(t₁)))
    #int_neurogenesis_2[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- param_bt[:wt].(times)))
    rawsim = simulate(valuemodel, Dict(:t₁ => Float64(700)))
    int_neurogenesis_wt[i] = sum(rawsim[:,2] .* 2 .* 0.95 .* (1 .- hillfun(t₁).(rawsim[:,1], times)))
end
metrics_neurogen[!,"value_noself"] = int_neurogenesis_1 ./ int_neurogenesis_wt
plot(xlab="Intervention Time (days)", ylab="increase in TAPs produced over wt", title="Life-long (700 days) neurogenesis", yticks=([100, 105, 110, 115], ["100\\%", "105\\%", "110\\%", "115\\%"]))
plot!(itimes, int_neurogenesis_1 ./ int_neurogenesis_wt .* 100, lab="'beta' model")
#plot!(itimes, int_neurogenesis_2 ./ int_neurogenesis_wt .* 100, lab="'value' model")
#plot!(itimes, int_neurogenesis_wt .* step, lab="wt", lc=:gray)

Big comparison

names(metrics_neurogen)

8-element Vector{String}:
 "itime"
 "type"
 "beta_fullself"
 "beta_partself"
 "beta_noself"
 "value_fullself"
 "value_partself"
 "value_noself"

ticks = 0.9:0.05:1.07
p = plot(legend=:bottomright, 
         ylim=(0.9, 1.07), 
         xlim=(0, 700),
         size=(300, 200),
         ytick=(ticks, ["$(round(Int64, 100*x))\\%" for x in ticks]),
         xlab="Time (days)", ylab="TAPs produced over wt")
Plots.abline!(0, 1, lab="", lc=:black)
#@df metrics_neurogen plot!(:itime, :beta_selfinteraction, lab="`beta` model with interaction", lw=2, lc=1)
#@df metrics_neurogen plot!(:itime, :beta_selfinvariant, lab="`beta` model without interaction", lw=2, lc=2)
@df metrics_neurogen plot!(:itime, :value_partself, lab="IFN dependent self-renewal", lw=2, lc=1)
@df metrics_neurogen plot!(:itime, :value_noself, lab="IFN independent self-renewal", lw=2, lc=2)
savefig("plot/life-long-pop-self.svg")

neurogen = @linq metrics_neurogen |>
    select(:t = :itime, :neurogen_dependent = :value_partself, :neurogen_independent = :value_noself)
loss = @linq metrics_loss |> 
    select(:t = :itime, :loss_dependent = :value_partself, :loss_independent = :value_noself)
innerjoin(neurogen, loss, on="t") |> CSV.write("metrics.csv")

"metrics.csv"

p = plot(legend=:topleft, xlab="Time (days)", ylab="Stem cells lost at 660 days",
         size=(300, 300),
         xlim=(0, 700), ylim=(0, 10))
#@df metrics_loss plot!(:itime, :beta_selfinteraction, lab="`beta` model with interaction", lw=2)
#@df metrics_loss plot!(:itime, :beta_selfinvariant, lab="`beta` model without interaction", lw=2)
#@df metrics_loss plot!(:itime, :value_selfinteraction, lab="Independent self-renewal", lw=2)
#@df metrics_loss plot!(:itime, :value_selfinvariant, lab="Combined self-renewal", lw=2)
@df metrics_loss plot!(:itime, :value_partself, lab="IFN dependent self-renewal", lw=2, lc=1)
@df metrics_loss plot!(:itime, :value_noself, lab="IFN independent self-renewal", lw=2, lc=2)
savefig("plot/loss-pop-self.svg")

Blub

struct ChangedBetaSelf <: ODEModel
    tspan::Tuple{Float64, Float64}
    times::Vector{Float64}
    values::Vector{Int64}
    fixed::Vector{Pair{Int64, Float64}}
    ChangedBetaSelf(tspan, times, values, fixed::Vector{Pair{I, N}}) where {I <: Integer, N <: Number} = new(tspan, times, values, fixed)
end
function ratefun(model::ChangedBetaSelf)
    function(du, u, p, t)
        _, _, βr, kab, nb_wt, nb_ko, pₛ, t₁ = p
        Q = u[1]; A = u[2]; r = u[3]
        b = 1 - hill(kab, t < t₁ ? nb_wt : nb_ko, Q+A)
        du[1] = dQ = -r * Q + 2* b * pₛ * A
        du[2] = dA = r * Q - pₛ * A
        if t < t₁
            du[3] = dr = -βr * r
        else
            du[3] = dr = 0
        end
    end
end
function initial(t::ChangedBetaSelf, x::AbstractVector)
    nsc₀, r₀, _, kab, nb_wt, _, pₛ, _ = x
    b₀ = hill(kab, nb_wt, 0)
    ratio = sqrt(((pₛ - r₀)/(2*r₀))^2 + (2*b₀*pₛ) / r₀)
    vcat(nsc₀ .* [1-1/(ratio+1), 1/(ratio+1)], [r₀])
end
link(t::ChangedBetaSelf, x::AbstractArray) = hcat(x[:,1] .+ x[:,2], x[:,2] ./ (x[:,1] .+ x[:,2]))
parameter_names(t::Type{ChangedBetaSelf}) = [:nsc₀, :r₀, :βr, :kab, :nb_wt, :nb_ko, :pₛ, :t₁]
bounds(t::ChangedBetaSelf) = [(100.0, 10000.0), (0.0, 1.0), (0.0, 0.1), (0.0, 1000.0), (0.01, 0.1), (0.01, 0.1), (0.0, 1.0), (0.0, 700.0)]
output_names(t::ChangedBetaSelf) = ["qNSC", "aNSC", "r"]
link_names(t::ChangedBetaSelf) = ["Total NSC", "Fraction active NSC"]

link_names (generic function with 4 methods)